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Abstract 


In the hydraulics industry, an accurate numerical approximation of the equations 
governing sediment transport in coastal regions has recently become a major topic of 
interest. These equations comprise the shallow water equations governing the water 
flow with the addition of a bed transport equation. It is common practice in industry 
to simplify the equations governing sediment transport by assuming the water flow 
is in an equilibrium state and the bed has a negligible effect on the water flow. 
However, this approach is limited as only steady flow and a slow moving bed can be 
approximated accurately in this manner. Thus, an unsteady approach is required 
that is considerably more robust and approximates the full system simultaneously. 
Until recently, the classic Lax-Wendroff scheme has been widely used in industry 
to obtain a numerical solution to the equations, but not surprisingly the numerical 
results obtained suffered from spurious oscillations resulting in the numerical scheme 
becoming unstable for long time periods. In industry, different measures have been 
applied to the use of classic Lax-Wendroff scheme to try to eliminate the spurious 
oscillations such as flux-limiter methods and using a small Courant number. 
Unfortunately, the spurious oscillations could not be eliminated and overpowered 
the numerical results for long computational run times. In this thesis, a variety 
of numerical schemes are discussed including adapted versions of the Lax-Friedrichs 
scheme, classic Lax-Wendroff scheme, MacCormack scheme and Roe’s scheme. High 
resolution schemes are also derived that satisfy the TVD property so that no spurious 
oscillations will occur in the numerical results. Five different formulations are then 
derived, which are based on either a steady or unsteady approach, and are used 
with the most accurate scheme and compared with the classic Lax-Wendroff scheme 
to obtain a numerical approximation of the equations. The numerical results are 
compared to determine which approach and numerical scheme is the most accurate. 
Sediment transport in both one and two dimensions is considered and dimensional 


splitting schemes are also discussed in two dimensions. 
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Chapter 1 


Introduction 


1.1 Motivation 


In recent years, sand transport has become a major topic of interest in the hydraulics 
community. The understanding of how sand interacts in certain environments is 
crucial for both the environment and businesses. For example, sand transport 
influences how harbours are constructed since if too much sand enters a harbour, 
the walls may become severely damaged and the costs of dredging the harbour may 
become too expensive and impractical. Reservoirs can lose all storage capacity due 
to sediment build up. Cunge et al. [3] stated that by 1973, 33% of the US reservoirs 
built before 1935 had lost between 25% to 50% of their original capacity due to 


sediment build up. 


In this thesis, we discuss a variety of numerical techniques that can be used to 
obtain a numerical solution of the equations that govern sediment transport. For a 


one dimensional channel, these comprise the equation for conservation of mass, 


Oh O(uh) 
Ot Ox 


=0, (1.1) 


Figure 1.1: The shallow water domain. 


the equation for conservation of momentum, 


O(uh) . O[hu? + kgh?] 
=—ghB,, 1.2 
a Oe ’ = 
and the bed-updating equation, 
OB Oq 
OE Sar ia 0, (1:3) 


where € = ae € being the porosity of the bed, which is non-dimensional, with 
0 <e< 1. In this thesis, we use a value of 0.4 for the porosity. Here n(z, t) 
represents the surface elevation, h(x, t) is the total height above the bottom of the 
channel (m), B(ax,t) is the height of the riverbed (m), u(a,t) is the velocity in 
the x direction (m/s) and q(u,h) is the total (suspended and bedload) volumetric 
sediment transport rate in the x direction (m?/s), see Figure 1.1. As suggested by 
the notation, the sediment transport flux q(u,h) is not in general a direct function 
of B, which can cause difficulties in obtaining an accurate numerical approximation. 
In some cases, the sediment transport flux cannot be written analytically and 
is calculated by using a “black box” approach where the flux is deduced from 
experimental data. Unfortunately, the equation for conservation of momentum (1.2) 
is not valid when discontinuities appear in the bed thus, a modified version of the 


equation is derived for this special case in Section 1.3. 


Until recently, the system (1.1) to (1.3) has been approximated by using a steady 
approach, which was pioneered by Cunge et al. [3]. The approach assumes that 
the changes in the bed have a negligible effect on the water flow and decouples 
the system. The system is decoupled into a water flow approximation, which 
is iterated to an equilibrium state, followed by a bed update. The approach is 
used with the classic Lax-Wendroff scheme in industry, but the scheme suffers from 
dispersion resulting in spurious oscillations occurring in the numerical results. 
Various techniques, including flux-limiter methods, have been used to try to eliminate 
the spurious oscillations. Unfortunately, the spurious oscillations could not be 
eliminated and overpowered the numerical results for long computational run times, 
see Damgaard [4] and Damgaard & Chesher [5]. Thus, even with the various 
techniques, the classic Lax-Wendroff scheme could only be used for short 
computational run times otherwise, spurious oscillations overpowered the numerical 
results resulting in the scheme becoming unstable. In this thesis, we discuss a 
variety of numerical techniques for approximating the equations governing sediment 
transport so that an accurate solution with no numerical oscillations present may 
be obtained. 


In this chapter, we discuss the derivation of the equations governing sediment 
transport together with two alternative sediment transport flux formulae. In 
Chapter 2, a variety of numerical schemes including adaptations of the classic 
Lax-Friedrichs, classic Lax-Wendroff, MacCormack and Roe’s scheme are discussed 
and applied to the shallow water equations without bed movement. To determine 
which scheme is the most accurate, three test problems are used to compare the 
numerical results of the different schemes. In Chapter 3, five different formulations 
of the governing equations are derived, which can be used to obtain a numerical 
approximation of the equations governing sediment transport. One of the 
formulations is based on the steady approach and the other four on the unsteady 
approach, where the water flow and bed update are approximated simultaneously. 
The flux-limited version of Roe’s scheme and the classic Lax-Wendroff scheme are 
then used to obtain a numerical solution of the different formulations for two test 
problems. The different formulations and numerical schemes are then compared to 


see which is the most accurate combination. In Chapter 4 the two dimensional 


case is considered and a variety of two dimensional schemes are discussed and 
applied to the 2D shallow water equations without bed movement. A test problem 
is used to determine which of the schemes is the most accurate. In Chapter 5, 
we discuss how to obtain an accurate approximation of the equations governing 
sediment transport in two dimensions. ‘The different formulations used in one 
dimension are adapted to two dimensions and are numerically approximated using 
the classic 2D Lax-Wendroff scheme, a flux-limited 2D version of Roe’s scheme 
and a dimensional splitting scheme. A conical sand dune test problem is used to 
determine which formulation and numerical scheme produced the most accurate 
numerical results. In Chapter 6, we discuss the results found in this thesis and 


propose further work. 


1.2 Derivation of the Equations 


We first derive the equations that govern fluid flow with sediment transport. By 
considering a region between x, and x2 in the one dimensional channel illustrated 


in Figure 1.1, we can derive the equations governing sediment transport. 


1.2.1 Conservation of Mass 


In the region x; to #2, we can determine that 


i. volume of fluid ial _ Bes of change of total a 


the region x1 to x2 of fluid in the region x; to x9 


Now, the total volume of fluid in the region x, to x2 is 


x2 h+B x2 x2 
/ | aya = f (h+B-B) ax = f h dz 
XY B £1 ry 


and by differentiating with respect to t, we obtain 


Rate of change of total a d / = hd 
= — va 
of fluid in the region 21 to x2 dt Ja, 
Also, 
aa oneal " Gi. and bee on " =i 
fluid entering at x, fluid leaving at x2 
thus, 


rs volume of fluid into 


= (uh)s, — (uh)ao- 
the region x1 to X9 | hah) (uh) es 


Hence, we obtain the integral form of the equation for conservation of mass 


ce i 
— 72 — 0), 1.4 
at J, h dx+ [uh]. 0 (1.4) 


To obtain the differential form, we integrate (1.4) with respect to t over the interval 
[t1,t2], where to > th, 


x2 x2 tg 
/ Ne eda / By dex / fuhJ®? dt =0. 


ial ty 


Then, by assuming that h(x,t) and u(x,t) are differentiable functions and using 


to ae h 
h(ata) — h(z,ts) = f a dt and [uh\?? -/ O(uh) co 


ty 


& 72 Oh (uh) 


Since x1, Yq and 1, tg are arbitrary, we obtain the differential form of the equation 


we obtain 


for conservation of mass, 


oh O(uh) 
Ot Ox 


Note that the common practice of re-writing this equation as 


= 0. 


On ir O(u(n — B)) OB 


ot Ox Ot 
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for the case of a fixed bed is not valid for the moving bed case considered here since 
BF 0. 


1.2.2 Conservation of Momentum 


In the region x; to #2, we can determine that 


ke rate of change of 


| = [Force applied in x direction] . 
momentum in x direction 


Now, 
l f ch f a. pe 
hae rate of change o | = al | u dydx + (hu)», — (hu?) 
momentum in x direction dt Je, Jp 
a 
a / fi de Op ~Pax. 
dt J, 
Also 
h+B ©2 
[Pressure force on ends] = g / (y—(h+ B)) wy 
B ry 
1 h+B] v2 
= 9 [By = (h-+ By | 
B dy, 
1 “e 
= |—<gh? 
aan 
and 
[Pressure force from riverbed] = —g | h dy, 
Czi,2 


where C’g;2 is the path of the line integral and denotes the curve of the riverbed in 


the region x, to £2. By assuming that no discontinuities are present in the riverbed 


dB 


v2 
[Pressure force from riverbed] = —g / h a 


1 


dx, 


we obtain 


sg — a6 |” we 
[Force applied in x direction] = ~agh -—g hs, dx. 


Ly ry 


Hence, 


a 9 a lel i‘ dB 
a . uh dx + (hu")s, — (hu a =| an" g : hs, dx 


ry 1 


and by re-arranging, we obtain the integral form of the equation for conservation of 
momentum, 
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d ee ™ dB 
= = = — dz. i 
at J, uh d+ [hu + 5h? af he dx (1.5) 


pa 


To obtain the differential form of the equation for conservation of momentum, we use 
the same approach as we did for the equation for conservation of mass and assume 


that h(x,t) and u(z,t) are differentiable functions and obtain 


O(uh)  O [hu? + Zgh?] 
Ot Ox 


= —ghB,. 


1.2.3. Bed-Updating Equation 


In the region x; to #2, we can determine that 


ia flux of mass i _ las of change of total mass 


the region x1 to x2 in the region x; to x 


Now, the total volume of sediment in the region x, to 22 is 


x2 B x2 
: i dydx = / B dz 
Ly 0 pan 


and by differentiating with respect to t, we obtain 


Rate of change of total mass a 
/ B ae. 


in the region x; to x2 


In addition we have 


hee volume of Total volume of 


sediment entering at 71 


| = €q(u,h),, and | = €q(u, h)oxo, 


sediment leaving at x2 
where € = ;+- and € is the porosity of the bed material, see Cunge et al. [3]. Thus, 


is flux of mass into 


| = €(q(u, h)x, — g(u, A)a») - 


the region x1 to X9 
Hence, we obtain the integral form of the bed-updating equation, 


a. “fre x 
— B dx+&[q(u, h)],? =0. (1.6) 
aE x, ’ 
To obtain the differential form of the bed-updating equation, we use the same 
approach as we did for the equation for conservation of mass and assume that h(x, t) 


and u(z,t) are differentiable functions and obtain 


1.3. Modification for a Discontinuity in the 
Riverbed 


When a discontinuity is present in the riverbed, see Figure 1.2, the integral 
conservation laws (1.4) and (1.6) are valid for non complex discontinuities 
(Needham & Hey [29] and Zanré & Needham [46]). However, the integral form 
of the equation for conservation of momentum (1.5) becomes invalid along the face 


of the discontinuity due to the term on the right side of (1.5) being derived with 


Figure 1.2: A discontinuity present in the riverbed. 


the assumption that no discontinuities are present in the riverbed. This term arises 


from the hydrostatic pressure force from the riverbed 


where C'g1,2 is the path of the line integral and denotes the curve of the riverbed in 


the region x; to 2. Away from the discontinuity, the line integral along C'g1,2 can 


-o | h dy =o [ ee di 
Cz1,2 Ly dx 


and equation (1.5) can be used. However, along the face of the discontinuity, we 


be determined 


cannot evaluate this integral as h(a,t) is undefined. Zanré & Needham [46] noted 


that when a discontinuity appears in the riverbed, we can use 
h=d'(y), min(Br, Br) < y < max(B_, Br), (1.7) 


where 
(ha —hy)(y— Br) 
Br—- By, 


d'(y) =hyt+ 


to determine the value of h(x,t) along the face of the discontinuity. By integrating 
(1.8), we obtain 
y(hr — hr)(y — 2Br) 


d(y) = hry + Br - Bp) 


(1.9) 


Thus, away from the discontinuity, the integral conservation laws (1.4), (1.5) and 


(1.6) can be used but at the discontinuity, (1.5) must be replaced with 


1 ey 
— uh dx + [h’ + an’ = -o | h dy 

2 x1 Cpi,2 
where (1.7), (1.8) and (1.9) provide a definition of depth along the face of the 


discontinuity 


1.4 Sediment Transport Flux Formulae 


The total load sediment transport flux includes suspended and bed-load sediment 
transport. Bed-load sediment transport includes the effects of grains of sand being 
transported on the surface of the bed by friction and gravity (for slopes). Suspended 
sediment transport includes the effects of grains of sand being picked up by the 
water flow and transported above the bed. For slow water flow, bed-load sediment 
transport is more dominant as not much sediment is carried by the water flow but 
as the water flow increases, suspended sediment transport becomes more dominant 
and sediment can be transported several meters above the bed especially if the grain 


size is small (see Soulsby [386] for more details). 


Numerous analytical sediment transport flux formulae have been derived that 
include both suspended and bed-load sediment transport and the choice of which 
to use is usually determined by the situation being modelled. In this thesis, we 


consider two of the more well known sediment transport flux formulae: 


e Grass [14] discussed one of the most basic sediment transport fluxes which 
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follows a simple power law, 
q(u) = Aulul”*. (1.10) 


Here, A is a dimensional constant (s?/m), that encompasses the effects of grain 
size and kinematic viscosity and is usually determined from experimental data 
with m being chosen so that 1 <m < 4. Unless wu is not allowed to change sign, 
(1.10) cannot be differentiated with respect to u. However, by considering odd 
integer values of m only, (1.10) can be differentiated and is valid for all values 


of u. For example, if we let m = 3 then 
g(u) = Aulu?| = Au, (1.11) 


which can now be differentiated with respect to u and is valid for all values of 


U. 


e Van Rijn [42, 43] derived a more complex sediment transport flux, 


(1.12) 


Au (|u| — ter)?" if Jul > Ue 
q(u, h) = . 
0 otherwise 


where U,, is the threshold current speed which is calculated from 


2h 
0.19(ds9)°" logio (=) if 100 < ds < 500m 
50 


Ucr = ’ 


2h 
8.5(ds0)°* logy (=) if 500 < dso < 2000,m 


dso (0.005 (42)"* + 0.012D,°) 
A= 
(gdso(s* — 1))'* 


(1.13) 


and the particle diameter is 


% 1/3 
D, = dso (Se — 0) . 
ay 


Here, dso is the median grain diameter in m, s* = ri is the specific gravity, ps 


tL 


is the mineral density of the sediment in kg/m* (2650kg/m? for quartz), p is 
the density of the water in kg/m? and 7¥ is the kinematic viscosity of the water 
in m?/s. Values of dso, ps and y can be obtained from Table 1.1, Table 1.2 
and Table 1.3 respectively. When calculating the threshold current speed, dso 


must be in m and not in jum. 


The sediment transport flux (1.12) is a simplified version of the full method of 
van Rijn [42] and is only valid for 0.5 < |u(z,t)| < 2.5m/s and 1 < h(z,t) < 
20m. Van Rijn derived (1.12) for fresh water, whose salinity is 0 parts per 
thousand (ppt), at a temperature of 15°C but the method can be used with 
a different salinity (i.e. salt water whose salinity is 35 ppt) and temperature 


with a minimal loss of accuracy. 


Notice that the sediment transport flux of van Rijn is considerably harder to 


differentiate than the sediment transport flux discussed by Grass. A variety of 


other sediment transport fomulae can be found in Soulsby [36] and van Rijn [43]. 


The sediment transport flux of van Rijn can be re-written in the form of (1.10) 


by assuming that the threshold wave speed is u,; = 0 and setting m = 3.4. Thus, by 


assuming that h is constant, we can use (1.13) to obtain a realistic approximation 


of A for (1.10). For example, if the situation being modelled is for seawater, whose 


salinity is 35ppt, with a temperature of 10°C, the following values are obtained 


— 


Ps = 2650kg/m?, which is the mineral density of quartz. 

p = 1027kg/m*, which is the density. 

s= S = 2.580331061, which is the specific gravity. 

v = 1.357 x 10-°m?/s, which is the kinematic viscosity. 

dso = 0.2 x 10-*m, which is the grain diameter of fine sand. 


h & 10m, which is the total height above the bottom of the channel. 


i 


Hence, by using the above values, we obtain the realistic value 
A & 0.00680396635637, 


which can be used with the sediment transport flux discussed by Grass. 


Mud / Sand dsy in um 
Fine Clay 0 to 0.9765625 
Medium Clay 0.9765625 to 1.953125 
Coarse Clay 1.953125 to 3.90625 
Very Fine Silt 3.90625 to 7.8125 
Fine Silt 1.8129 to 15.625 
Medium Silt 15.625 to 31.25 
Coarse Silt 51.25 to 62.5 
Very Fine Sand 62.5 to 125 
Fine Sand 125 to 250 
Medium Sand 250 to 500 
Coarse Sand 500 to 1000 
Very Coarse Sand 1000 to 2000 


Table 1.1: The grain size in jum (taken from Soulsby [36]). 


1.5 Summary 


In this chapter, we have derived the equations governing sediment transport in one 
dimension. We have also discussed two different sediment transport formulae that 
can be used. In the next chapter, we discuss a variety of numerical schemes that can 


be used to numerically approximate the equations governing sediment transport. 
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Temperature °C 


Salinity 0 5 10 15 20 25 30 39 40 
0 999.8 1000 | 999.7 | 999.1 | 998.2 997 995.7 994 992.2 
5 1003.9 | 1004 | 1003.6 | 1003 1002 | 1000.8 | 999.4 | 997.7 | 995.9 
10 1008 | 1007.9 | 1007.5 | 1006.8 | 1005.8 | 1004.6 | 1003.1 | 1001.4 | 999.6 
15 1012 | 1011.9 | 1011.4 | 1010.6 | 1009.6 | 1008.3 | 1006.8 | 1005.1 | 1003.2 
20 1016 | 1015.8 | 1015.3 | 1014.4 | 1013.4 | 1012.1 | 1010.5 | 1008.8 | 1006.9 
25 1020 | 1019.8 | 1019.2 | 1018.3 | 1017.2 | 1015.8 | 1014.3 | 1012.5 | 1010.6 
30 1024.1 | 1023.7 | 1023.1 | 1022.1 | 1021 | 1019.6 | 1018 | 1016.2 | 1014.3 
30 1028.1 | 1027.7 | 1027 1026 | 1024.8 | 1023.3 | 1021.7 | 1019.9 | 1018 
40 1032.2 | 1031.6 | 1030.9 | 1029.8 | 1028.6 | 1027.1 | 1025.5 | 1023.7 | 1021.7 


Table 1.2: Density of water in kg/m* (taken from Ramsing & Gundersen [32]). 


Temperature °C 

Salinity 0 5 10 15 20 25 30 35 40 
0 1790.1 | 1517.5 | 1305.2 | 1140.2 | 1009.9 | 901.5 | 802.2 | 699.5 | 580.4 
5 1795.2 | 1524.5 | 1312.9 | 1147.8 | 1016.9 | 907.7 | 807.9 | 705.1 | 586.9 
10 1800.2 | 1531.4 | 1320.5 | 1155.3 | 1023.8 | 913.9 | 813.6 | 710.8 | 593.4 
15 1805.3 | 1538.3 | 1328.1 | 1162.8 | 1030.7 | 920 | 819.2 | 716.3 | 599.8 
20 1810.3 | 1545.2 | 1335.6 | 1170.2 | 1037.5 | 926.1 | 824.7 | 721.9 | 606.2 
25 1820.2 | 1558.7 | 1350.5 | 1184.8 | 1050.9 | 938.1 | 835.7 | 727.4 | 612.5 
30 1820.2 | 1558.7 | 1350.5 | 1184.8 | 1050.9 | 938.1 | 835.7 | 732.8 | 618.8 
35 1825.1 | 1565.4 | 1357.8 | 1192 | 1057.6 | 944.1 | 841.1 | 738.2 | 625 
40 1829.9 | 1572 | 1365.1 | 1199.2 | 1064.2 | 949.9 | 846.5 | 743.6 | 631.2 


Table 1.3: Kinematic viscosity of water in x10~°m?/s (taken from Ramsing & 
Gundersen [32]). 
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Chapter 2 


Numerical Schemes for Systems of 
Conservation Laws in One 


Dimension 


Now that we have derived the equations governing sediment transport, we discuss 
a variety of numerical schemes that can be used to numerically approximate the 
equations. We discuss an adapted version of the Lax-Friedrichs (LxF') scheme (see 
Garcia-Navarro et al. [9]), the staggered, non-staggered and central LxF and NT 
scheme (see Nessyahu & Tadmor [30] and Jiang et al. [22]), the Lax-Wendroff scheme 
(see Lax & Wendroff [24]), an adapted version of the MacCormack approach (see 
LeVeque & Yee [27], Yee [45] and Hudson [20]) and an adaptation of Roe’s Scheme 
(see Roe [33], Glaister [11] and Hubbard & Garcia-Navarro [18]). Flux limiter 


methods and slope limiter methods are employed to minimise numerical oscillations. 
The equations governing sediment transport can be written as a system of 
conservation laws with source term, i.e. 


Ow OF (w) 
Ot Ox 


=h. (2.1) 
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where F(w) is the flux-function and R is the source term. Before considering the 
full system, we illustrate the numerical techniques for systems of conservation laws 


using the example of the shallow water equations, 


ia) 


which are a system of conservation laws. Some of the numerical schemes discussed 


0 


: 2.2 
_ghB, (2.2) 


uh = 
hu? + sgh? 


in this chapter require the Jacobian matrix associated with the system (2.1), which 


for the shallow water equations is 


whose eigenvalues are 


My =u-—Vgh and do = ut >/gh, 


with corresponding eigenvectors 


2.1 Test Problems 


We use the following three test problems for the shallow water equations to illustrate 
the accuracy of the numerical schemes discussed in this chapter. For all three test 


problems, the riverbed is fixed. 
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2.1.1 Test Problem A: The Dam Break Problem 


For this test problem, the riverbed is of constant depth thus, the source term is no 
longer present, R = 0. The test problem consists of a 1D channel of length 1m with 
walls at either end. The initial velocity is 0 and a barrier is present at « = 0.5, 


which is removed at t = 0. The initial conditions consist of 
u(z,0)=0 and. Atz,0)= 


and are illustrated in Figure 2.1. 


h A 
DAM 


hy 


hr 


QY 


0 0.5 1 


Figure 2.1: Initial conditions of the dam break problem. 


Stoker [38] derived an analytical solution of the dam break problem which can 


be used to illustrate the accuracy of the numerical schemes. The exact solution is 


1 
ula,t) =) 3g C+ tVahn) 1) ify —tVohy Sas (ue) +3 
| 2 if (uz —c)t +5 << St+3 
0 if > St+ 4 


AW 


and 


hy ife<$—tV/ghr 
1 1 ° 
Oa (2vari — 5 (22 — 0) if §—t/ghy <2 < (ue—e)t+ 3 


g 
= 89 
2 ghr 


hr fas St+2 


where 


_ ghr 8S? = ghr 8S? 
us = 5 18 G ari and. = 5 a a 1]. 


The bore speed, S, is the positive root of 


Ug + 2co — 2\/ghy, = 0. 


The exact solution with hz = 1m and hz = 0.5m is illustrated in Figure 2.2 and 
Figure 2.3, where the bore speed is approximately S = 2.957918120187525. The 
exact solution is only valid until either the bore or the rarefraction wave hits the 
walls. For a more in depth analysis of the exact solution, see Stoker [38] and 
Glaister [11]. 


2.1.2 Test Problem B: Wave Propagation Test Problem 


LeVeque [26] discussed a wave propagation test problem with a pulse present in the 
riverbed and a small disturbance in the river. Since the riverbed is not constant, a 


source term is now present. The initial conditions are 


l+w—B(x) if01l<2<0.2 
1— Biz) otherwise 


ue,0)=0, he0)= 
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h(x,t) 


Oo 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 


Figure 2.2: 


O------ 0.02 0.04 —S— 0.06 —x— 0.08 0.1 


Figure 2.3: Exact solution of the dam break test problem for t = 0 to 0.1s (wu). 
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EF BX) 


h(x,0) + B(x) 


Figure 2.4: Initial conditions of the wave propagation test problem with w = 0.2 (h 
& B). 


and the bathymetry is 


rs 


Bt + (cos (107 (a — $)) +1) if |je-sl<a 
aa fice: otherwise 


which are illustrated in Figure 2.4 for w = 0.2. The value of w is taken as either 0.2 
or 0.01 and we follow LeVeque [26] and use a gravitational constant of g = 1. The 
disturbance in the river created by w splits into two waves propagating in opposite 
directions and for small w the characteristic speeds are approximately +./gh. For 


small w, an accurate numerical solution is harder to obtain. 
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—= B(x) —— h(x,0) + B(x) 


Figure 2.5: Initial conditions of the tidal wave propagation test problem (h & B). 


h(0,t) 
@ 
& 
a 


Oo 10800 21600 32400 43200 54000 64800 


Figure 2.6: Upstream boundary condition of the tidal wave propagation test problem 


(h). 
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2.1.3. Test Problem C: Tidal Wave Propagation Test Problem 


Bermtidez & Vazquez [1] discussed a tidal wave propagation test problem that 


consists of the initial conditions 
u(z,0)=0, A(x,0) = 60.5 — B(x) 


and bathymetry 


B(x) = —* +10 (1-sin(a(F+5))), 


where L is the length of the channel which is usually taken as L = 648,000m. The 


physical boundary condition, 


(0,1) = 64.5 + 4sin (4 (755 -1)) if t < 43, 200s 
60.5 if t > 43, 200s 


is used to simulate a tidal wave of 4m amplitude entering the region at the upstream 
boundary, see Figure 2.6. At the downstream boundary, the physical boundary 
condition u(L,t) = 0 is also required. The tidal wave reaches a full height of 8m at 
t = 21, 600s and since the wave propagates at approximately /gh, at t = 10, 800s 
the tidal wave should have only reached as far as x © 216,000m for L = 648, 000m. 


2.2 Conservative Numerical Schemes 


We can accurately numerically approximate a system of conservation laws by using 


a conservative numerical scheme with source term approximation, 


wt! — w" — 6(Ft_, — Ft_1)+sRf, (23) 


ae 1 
da 3 
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Figure 2.7: The one dimensional mesh. 


_ At : : 3 : Q 3 
where s = <,, Az and A? are the step sizes in space and time respectively, F* 41 1s 


the numerical flux, R# is the source term approximation and 


1 "itd 
w) & — wat") da 
ie 
Lv. 4 
= 
2. 
is the numerical approximation, see Figure 2.7. Here, the points x = x_1 and 
2 
© = £,,1 are the spatial boundaries and we require numerical boundary conditions 
2 


at these points. The spatial step size, Az, is fixed and we use a variable time step, 


vAL 
max;(|Az|)’ 


Ai= 


where A; are the eigenvalues of the Jacobian matrix and v is the required CFL 
number (Courant, Friedrichs and Lewey). Unless stated, all schemes discussed in 


this chapter are stable for v < 1. 


To ensure the error of a numerical scheme does not grow, the variables are non- 
dimensionalised so that the spatial and time step-sizes are less than one, i.e. Ar < 1 


and At < 1. For the shallow water equations, we non-dimensionalise the variables 
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by using 


where 


L 
L=|x;—2%| and T= 4/— 
g 


denote the non-dimensional coefficients and L is the length of the domain. 


2.2.1 Source Term Approximation 


In gas dynamics, source terms can become stiff being extremely difficult to 
approximate accurately, see LeVeque & Yee [27] for more details. For the shallow 
water equations, the source term becomes increasingly difficult to approximate as 
the variation in the riverbed becomes more pronounced, see Hudson [20] for more 


details. In general, the source term can be approximated in two ways: 


e a pointwise approach, where the source term approximation is calculated at 
the nodal points, i.e. 
R} = ArcR(w7). 


e an upwind characteristic based approach, where the source term is 
approximated in a more physical way by averaging the source term, 
_ Ag 


Bia > (Rhy +REy) 


and obtaining upwind approximations, R?,, of the source term. 


dle 


In general, the pointwise approach is considerably less accurate than the upwind 


approach due to the pointwise approach not satisfying the C-property, which is 


discussed in the next section. 
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Figure 2.8: The quiescent flow case: u = 0 and h= D— B. 


2.2.2 C-Property 


Bermudez & Vazquez [1] and Vaézquez-Cendén [44] discussed an approach for 
approximating source terms which is designed for quasi-steady and steady flow. 


Consider the shallow water equations for the quiescent flow case, 
u(z,t)=0 and A(z,t)=D-— Bi(az,t) V(z,t), 


see Figure 2.8. For this stationary case w; = 0 and thus the flux function and source 
term balance: 
Fy => fh. 


Therefore, an accurate numerical scheme would also balance the numerical flux with 


the source term approximation, 
* * _ -R* 
Fiya eee =i, 


Hence, we derive numerical schemes that satisfy the C-property (conservation 
property) when applied to the quiescent flow case. If the source term approximation 


balances with the numerical fluxes, then the numerical scheme satisfies: 


e the approximate C-property, if the numerical scheme is accurate to the order 


O(Ax?) when applied to the quiescent flow case; 
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e the exact C-property, if the numerical scheme is exact when applied to the 


quiescent flow case. 


If a numerical scheme does not satisfy the C-property (exact or approximate) then 


spurious waves may occur in the numerical results. 


2.3 High-Resolution Schemes 


Unfortunately, second order numerical schemes suffer from dispersion resulting in 
spurious oscillations appearing in the numerical results, especially around 
discontinuities. This is due to second order numerical schemes not satisfying the 
Total Variational Diminishing (TVD) property, which is discussed in the next 
section. Thus, we construct a numerical scheme that satisfies the TVD property 
by using a second order accurate numerical scheme on smooth solutions and adding 
diffusion to the numerical scheme near discontinuities. We call such numerical 
schemes high-resolution schemes, which are at least second order accurate on smooth 


solutions and minimise the spurious oscillations present near discontinuities. 


2.3.1 Total Variational Diminishing 


Consider the scalar conservation law 
which can be numerically approximated by using the conservative numerical scheme 


wht! =wP—s (fia - fia). (2.5) 
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The exact solution of a scalar conservation law (2.4) has Total Variation, 


TV(w) = f fwsl ae, 


which does not increase and only decreases across shocks, see Lax [23] for more 
details. For scalar conservation laws (2.4), this property can be used to eliminate 
spurious oscillations present in a numerical scheme. Unfortunately, this property 
does not hold for systems of conservation laws but the design criteria based on 
this property can still be used to minimise the spurious oscillations present in a 
numerical scheme when applied to systems. Harten [15] deduced that a numerical 


scheme should also satisfy the Total Variation Diminishing (TVD) property, 


TV(w"t") < TV(w"), 
where the numerical Total Variation is defined by 


TV(w"tt) = So writ — w?t], 
a 


For scalar conservation laws, if a numerical scheme satisfies the TVD property, then 
no spurious oscillations will occur in the numerical results. Harten [15] discussed an 


approach that determines if a numerical scheme is TVD by re-writing (2.5) as 


BE a yy? =f. Di, (we) - wr?) - Cia (w? - tii) : 


w i 


and determining if the inequalities 


are satisfied Vi. If the inequalities are satisfied, the numerical scheme satisfies the 


TVD property. 


ot 


Flux-limiter &(0) 


Minmod max(0, min(1,6)) 


Roe’s Superbee | max(0, min(26, 1), min(0, 2)) 


van Leer Le 
1+ |6| 
67 +6 

Alb 
van ada Te 


Table 2.1: Some flux-limiters 


2.3.2. Flux-Limiter Methods 


One approach we can use to construct a high-resolution scheme is to use a flux-limiter 
method as discussed by Sweby [40] and LeVeque [25]. For the scalar conservation 
law (2.4), we used the conservative numerical scheme (2.5). Flux-limiter methods 


construct a numerical flux of the form 


fie? = fF9 + 05, (685 - 9), (2.6) 
where nee is a second order numerical flux, raw is a first order numerical flux and 
Or is a limiter. The limiter is chosen such that if the data is smooth, then the 
limiter is 1, which reverts the numerical flux to a second order approximation but if 
the data is near a discontinuity, the limiter is 0, which reverts the numerical flux to 
a first order approximation. van Leer [41] and Roe [34] discussed an approach we 
can use to measure the smoothness of the data by looking at the ratio of consecutive 
gradients, 7 : 

C4 = ie ae where J =i—sgn Ca ; 


Here, is the wave speed and if Ona is close to 1, then the data is smooth, but 


if #1 is near 0, then a kink is present in the data. Thus, we take ®”, to be a 
2 2 
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Figure 2.9: A piecewise constant and linear representation. 


function of 0” ,, i.e. 
i+5 
ae =o (07, 1) ’ 
where ® is a given function. Sweby [40] discussed a variety of limiters that guarantee 


second order accuracy whilst still satisfying the TVD property, which are listed in 
Table 2.1. 


2.3.3 Slope-Limiter Methods 


We can also use a geometric approach to obtain a high-resolution scheme. 
LeVeque [25] discussed slope-limiter methods in detail and noted that to obtain 


a slope-limiter method, we must construct a numerical scheme as follows: 


1. We use the numerical data w} and construct a piecewise linear function 


w"(x,t") = we +o7(@—2;) where t1<t< 241 


and 0? = au is the slope of the 7th cell, see Figure 2.9, which is based on the 


numerical data. 
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2. The exact solution of a generalised Riemann problem, w”(x,t”*'), which is 


related to (2.4) is then obtained using the piecewise linear function. 


3. The exact solution is then averaged, 


n+l l Mid ~n n+l 
lea er 0 ae 


ns 
ame?) 


n+1 


to obtain the new data w;’"~ and the process is repeated. 


Consider the following example, as discussed by LeVeque [25], for the scalar linear 
advection equation, 


wt + aw, = 0, 


whose exact solution is 
w(a,t) = w(x — at, 0). 


By constructing the piecewise linear function described in step 1 from the data, we 


obtain the exact solution 
w'(a, f°") = 8" (a — at, t”), 
and by integrating the exact solution (for a > 0) as described in step 3, we obtain 
wp = wl — (wf — why) — 5aAt(1 —v)(o? — of). (2.7) 


Now, we must choose the slopes, o/’, such that the scheme (2.7) is second order 
accurate and satisfies the TVD property. If we choose o? = 0, we obtain the first 
order upwind scheme and by choosing a? = x (wry — w!'), we obtain the classic 
Lax-Wendroff scheme. If the upwind slopes are used, the scheme is not second order 
accurate and if the Lax-Wendroff slopes are used, then overshoots may occur in 
the piecewise linear representation, see Figure 2.9, which results in an increase in 
Total Variation. Thus, we can view these as being a poor choice of slopes. A better 


choice of slopes, which makes the scheme second order accurate and satisfy the TVD 
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property, is to use the minmod limiter, 


oss 
of = —minmod(w?,, — w?,w? — wP,), 


: Ax 


where 


minmod(a, b) = 5 (sen(a) + sgn(b))min(|a], |b]). 


It is also interesting to note that by setting 


a n n n 
= Tg wits ~ Wi )Pi 2, 


the scheme reverts back to a flux-limited scheme, where ©” , is a flux-limiter, which 
2 


was discussed in the previous section. Thus, the minmod limiter can also be used 


as a slope-limiter. 


2.4 Boundary Conditions 


Boundary conditions are required at the upstream and downstream boundaries. For: 


e Test Problem A, walls are present at the upstream and downstream boundaries 


and thus we need to reflect the velocity 


n _ n no n 
ue, = ue, and uly, = up 44 


—2 


and for the water depth, we use the simple boundary conditions 


where 7 = 1,2. 


e Test Problem B, we require transmissive boundary conditions due to the 
disturbance in the river splitting into two waves travelling in opposite 


directions, which we discuss next. 
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e Test Problem C, we use the physical boundary conditions 
Ae— 0,0") and: ap. =0, 


where 7 = 0,1,2 and the numerical boundary conditions 


n _ £0 n _.n 
TH = hy and ut; = up, 


where 7 = 1, 2. 


2.4.1 'Transmissive Boundary Conditions 


In some test problems, waves pass through the boundaries and, if the incorrect 
boundary conditions are chosen, then the wave may be reflected back into the domain 
resulting in an inaccurate numerical solution. Transmissive boundary conditions 
allow waves to leave the region without being reflected back into the domain. One 
type of transmissive boundary conditions can be obtained by using the k-Riemann 


invariants (see Godlewski & Raviart [13]), which are obtained by solving 


where e; are the eigenvectors, r denotes the k-Riemann invariants and k is the k“” 


component of the system. For the shallow water equations we have 


aa t (u- v ah) sai =n 


and 


a t (ut von) aa = 


and by solving the two equations, we obtain 


rp =u-—2/gh and rg=ut+2/gh. 
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Hence, by using the k-Riemann invariants, we obtain the following approximations 
of u and h, 


1 1 
oS re) and h= 76, (721) 


We can use these approximations of u and h at the boundaries to allow waves to 
pass through without being reflected back into the domain. For subcritical flow, i.e. 
Vgh > |u|, r; and rz represent the left and right moving waves respectively. Thus, 
at the upstream boundary, we want to eliminate all the right moving waves but let 


the left moving waves pass through the boundary, which may be obtained by setting 
ry =u —2/f/ght and roe=ugt+2/ghj, 


where uo and hj denote the initial values of the velocity and height of the river at the 
upstream boundary respectively. Similarly, at the downstream boundary, we want 
to eliminate all the left moving waves but let the right moving waves pass through 


the boundary, which may be obtained by setting 
Ty = uz —2./gh, and z= u7_4+24/ ght; 


where u; and h7 denote the initial values of the velocity and height of the water 


above the bottom of the channel at the downstream boundary respectively. 


2.5 LxF and NT Scheme 


2.5.1 Adapted LxF Scheme 


One of the most basic numerical schemes we can use to approximate (2.1) is the 


Lax-Friedrichs (LxF) scheme as discussed by Garcia-Navarro et al. [9], 


i 
wotl — Ow? + —_— = 


5 (Fni- FL) +R}, (28) 


(wey + wy" 1) 
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where 0 < 6 < 1. The scheme is unstable for a value of 9 = 1 and as 6 —> 0, 
the scheme becomes more stable but unfortunately more diffusive. The value of 
§ = 0.1 is most commonly used and @ = 0 is the classic Lax-Friedrichs scheme. 
Unfortunately, the Lax-Friedrichs scheme suffers badly from diffusion and a 
“staircase” effect, where kinks are present in the numerical results, can 
sometimes occur in the numerical results, especially for 6 = 0. However a numerical 


approximation of the Jacobian is not required. 


We use a centralised pointwise approach to approximate the source term, R’, 


and for the shallow water equations this takes the form 


Unfortunately, we cannot derive a source term approximation that makes the Lax- 


Friedrichs scheme satisfy the C-property. For the quiescent flow case, 
1-0 
ne = ont +258 (np. + ny) 
and since h © D— B and h?'t! = h”, we obtain 


(0 — 1) (B41 — 2B; + Bis) = 0, 


which can only be satisfied if either 6 = 1 or the riverbed is linear. Taking 6 = 1 
renders the numerical scheme unstable and the source term is zero when the bed is 


linear. 


2.5.2 Staggered and Non-Staggered LxF Scheme 


Nessyahu & Tadmor [30] and Jiang et al. [22] discussed a first order central scheme 
called the staggered LxF scheme, 


n 1 nm n n n n 
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Figure 2.10: Staggered grid. 


which is a modified version of the classic Lax-Friedrichs scheme (LxF) and is derived 
on a staggered grid, see Figure 2.10. The scheme has been adapted to include a 
pointwise source term approximation, which for the shallow water equations takes 
the form 
R= : 
—_ | —S (hia + ht)(Be — Be) 


We can obtain a non-staggered version of this scheme by averaging the two 


neighbouring staggered cells, 


1 
a rr n+1 n+1 
W; = 5 (watt + w’ :) 


| 
z 
i 
oe 
bo 
E 
oe 
2 


s s 
41) =) (Fey - F? 1) gx 2 (Re. zi R! :) . (2.10) 


By averaging the staggered version of the scheme, we slightly reduce the accuracy 
but obtain numerical results at the desired points. As with the adapted LxF scheme, 
both the staggered and non-staggered version of the LxF scheme do not satisfy the 
C-property and are stable for v < 0.5. Notice that by using the staggered or non- 
staggered LxF scheme, the Courant number has halved compared to the adapted 


LxF scheme. 
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2.5.3 Central, Staggered and Non-Staggered NT Scheme 


Nessyahu & Tadmor [30] extended the staggered LxF scheme to second order and 
derived a high-resolution non-oscillatory central difference scheme called the 
staggered NT scheme, 


wtl = = (w2, +-w?)—= ((wa)2.4 — (we)?) —s (any = F *?) +eR™'2, (2.11) 
its ) i+1 a 8 i+1 a i+1 a its ’ 


where 
wi? = wh +5 (R} — (F)f). 


a 7 


The staggered NT scheme (2.11) has been adapted to approximate a system of 
conservation laws with source term present (2.1). To ensure that the staggered NT 


scheme is non-oscillatory, we use either the limiter 
* =MM(ad ! A yas A 212 
(w.); = SAW ek 5 ( Wap ay w;-1) pOEWeed ) 9 ( . ) 


where 0 < a < 4, or the UNO limiter [16] 


1 1 
(w..); = MM (a. - 53MM (A? wisi, A’w;) , Aw,_1 + 3MM (A?w;, A’wi-1) } , 


a) 
(2.13) 
where Aw;,1 = wh, —w?, A?w; = wt, — 2w? + w?,. Here, MM denotes the 


minmod limiter 


| dmin(la|, |, |el) if d = sgn(a) = sgn(b) = sgn(c), 
minmod(a, b,c) = 
0 otherwise. 
The value of (F,)? can be determined from either the UNO limiter, (2.12) or from 
(F,)? = A?(w,)?. Notice that setting the numerical approximations (w,)? = 


(F.,)!' = 0 reverts the scheme to the first order LxF scheme. 


We use a centralised pointwise approach to approximate the source term, R?, 


36 


and for the shallow water equations this takes the form 


0 
—gh? (Bz)? 


a 


A non-staggered version of the NT scheme can be obtained with a slight loss 
of accuracy by averaging the neighbouring staggered cells with a piecewise-linear 


interpolant to obtain 


1 i 
nt+1 __ n+1 n+1 i = 
Ww; = 9 (wr + wit!) = (wa)? _— (wa)? 1) ; 


8 
where 
(Wo)i,4 = MM (AwfT, Aw?'*") and Age = Wee - wer. 
Hence, 
wytt = 2 (wis + wi + wh) ((we)Baa — (walls) —g ((we)Pa — (wala) 
4 16 8 2 2 
-3 Gx = Pt) +5 a + RY?) , (2.14) 


which is the non-staggered NT scheme. Both the staggered and non-staggered 
version of the NT scheme are stable for vy < 0.5 and as with the first order LxF 
schemes, neither satisfy the C-property. 


Nessyahu & Tadmor [30] also extended the classic central Lax-Friedrichs scheme, 


(2.8) with 0 = 0, to second order. The high-resolution central NT scheme is 


n i n n l n n 5 nts n+5 n+s 
Ww; a 2 (wr + wis) | ((We) fe a (w:)i4) ~ 95 (aa - i) +sR; ’, 
(2.15) 
where 
w, ? = wy += (RP -(F,)f) 
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and for the shallow water equations the source term approximations are 


0 
Rr =| 


n+4 0 
: 7 and R,-? = g 
—ghi (Bz); 


nti nti n+i nti 
~ hin 7 hi (Big? — By’) 


The high-resolution central NT scheme is stable for vy < 1 and does not satisfy the 
C-property. If we set the numerical approximations (w,)? = (F,)? = 0, then the 


central NT scheme reverts back to the classic Lax-Friedrichs scheme. 


2.5.4 Numerical Results of the First Order LxF Schemes 


To compare the different versions of the LxF schemes, we use Test Problem A, which 
is the dam break test problem. The dam break test problem is used with hy = 1m 
and hy = 0.5m. The central and adapted LxF schemes are used with v = 0.8 and 
the staggered and non-staggered LxF schemes are used with vy = 0.4. By using 
the different first order LxF schemes with Ax = 0.01m, we obtain the numerical 
results in Figure 2.11 and Figure 2.12. The dam break test problem also has an exact 
solution which is illustrated to show the accuracy of the different LxF schemes. Here, 
we can see that the staggered LxF scheme produced the most accurate numerical 
results but the scheme was still very dissipative. The non-staggered LxF scheme 
produced similar numerical results to the staggered LxF scheme but they were less 
accurate due to the scheme being more diffusive. The adapted LxF scheme with 6 = 
0.1 was more diffusive than both the staggered and non-staggered LxF schemes. The 
central LxF scheme, which is the classic Lax-Friedrichs scheme, produced completely 


inaccurate numerical results due to the “staircase” effect. 


Test Problem A does not have a source term present thus, we use Test Problem 
B, which is the wave propagation test problem to determine which first order LxF 
scheme is the most accurate when a source term is present. For this test problem, 
we use w = 0.2 and Ax = 0.01. The central and adapted LxF schemes are used with 
vy = 0.8 and the staggered and non-staggered LxF schemes are used with vy = 0.4. A 


flux-limited second order version of Roe’s scheme, which is discussed in Section 2.7, 
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x 


Figure 2.11: Numerical results of the different first order LxF schemes for Test 
Problem A at t = 0.1s (h). 


u(x,0.1s) 


Figure 2.12: Numerical results of the different first order LxF schemes for Test 
Problem A at t = 0.15 (u). 
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on a fine mesh, Ax = 0.0001, is also used as a reference solution. From Figure 2.13 
and Figure 2.14, we can see that none of the different first order LxF schemes have 
produced accurate results compared to the fine mesh numerical results. All of the 
numerical schemes suffered badly from diffusion and have produced oscillations over 
the pulse in the riverbed. This is due to none of the numerical schemes satisfying 


the C-property and thus oscillations have occurred in the numerical results. 


2.5.5 Numerical Results of the High-Resolution NT Schemes 


To compare the different versions of the high-resolution NT schemes, we use Test 
Problem B, which is the wave propagation test problem. For this test problem, we 
use w = 0.2 and Ax = 0.01. The central and adapted NT schemes are used with 
v = 0.8 and the staggered and non-staggered NT schemes are used with v = 0.4. A 
flux-limited second order version of Roe’s scheme, which is discussed in Section 2.7, 
on a fine mesh, Ax = 0.0001, is also used as a reference solution. From Figure 2.15 
and Figure 2.16, we can see that as with the first order LxF schemes, none of the 
different high-resolution NT schemes have produced accurate results compared to 
the fine mesh numerical results. All of the numerical schemes suffered badly from 
oscillations and the central NT scheme has again produced the “staircase” effect. 
This is due to none of the numerical schemes satisfying the C-property and thus, 


oscillations have occurred in the numerical results. 


Hence, when a source term is present, the LxF and NT schemes produce 
inaccurate numerical results due to spurious oscillations being present. These 
spurious oscillations can be reduced by making the LxF and NT schemes satisfy 
the C-property but unfortunately, no source term approximation can be derived 


that makes the schemes satisfy the C-property. 
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Surface elevation 


Fine Mesh ------ Adapted LxF ———__ Central LxF —-G— Staggered LxF —-x— Non-Staggered LxF 


Figure 2.13: Numerical results of the different first order LxF schemes for Test 
Problem B with w = 0.2 at ¢=0.7 (h+ B). 


u(x,0.7) 


Fine Mesh -- --- - Adapted LxF —— Central LxF —@— Staggered LxF —+<— Non-Staggered LxF 


Figure 2.14: Numerical results of the different first order LxF schemes for Test 
Problem B with w = 0.2 at t = 0.7 (u). 
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Surface elevation 


Fine Mesh ------ Central NT —_ Staggered NT —-C— Non-Staggered NT 


Figure 2.15: Numerical results of the different high-resolution NT schemes for Test 
Problem B with w = 0.2 at t =0.7 (h+B). 


u(x,0.7) 


Fine Mesh - - - - -- Central NT ——- Staggered NT —-G— Non-Staggered NT 


Figure 2.16: Numerical results of the different high-resolution NT schemes for Test 
Problem B with w = 0.2 at t = 0.7 (u). 
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2.6 Lax-Wendroff & MacCormack Scheme 


2.6.1 Lax-Wendroff Scheme 


The Lax-Wendroff scheme [24] is derived from a Taylor’s series expansion, 


Ar 
wit! we wi + At(w,)? + a (wee)? + O(At?), 


where the partial differential equation is used to replace the time derivatives with 


spatial derivatives and leads to the Lax-Wendroff numerical flux-function 
n n 8 n nm va) 
(Fi. a F?) —~ 5fhigd (Fas = F’) : (2.16) 


The Lax-Wendroff scheme requires a numerical approximation of the Jacobian 


matrix, which can be approximated by averaging the neighbouring cells, 


Anp= a (TENE) (2.17) 


For the shallow water equations, we can use a pointwise approach, 


0 
| na ae Bey UB Be a) 


4 


(2.18) 


to approximate the source term but this is a crude approximation as it does not 
satisfy the C-property and reduces the accuracy of the numerical scheme as the 
scheme is no longer strictly second order accurate. Another more accurate approach, 
which maintains second order accuracy is to include the source term in the Taylor 


series expansion to obtain 


R* = AiR + —— (R, — (AR),). 
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Thus, we obtain 


il h . 7 . 
= 9 ((i-sAty) Bly + (T4582 y)RLy) (2.19) 


where for the shallow water equations, 


0 
R” ‘ro 
m4 | —5 (hiya + AP)(BE — BY) 


Here, the R; term has been omitted resulting in a slight loss of accuracy but the 
Lax-Wendroff scheme with (2.19) satisfies the C-property. Alternatively, we can 
approximate this term by using the chain rule R; = Ryw;. Thus, we obtain a 


semi-implicit version of the Lax-Wendroff scheme 
r = 
wrtlowt—s k 2 > (Rw)? ( es ee R;) (2.20) 
2 +5 13 


where for the shallow water equations 


0 0 | 
=3( Big = Bi 4) 0 


Unfortunately, the semi-implicit approach cannot be used for all systems of 
conservation laws as the derivative of R with respect to w can sometimes be difficult 


to obtain. 


2.6.2. Flux-Limited Lax-Wendroff Scheme 


The Lax-Wendroff scheme suffers from dispersion resulting in spurious oscillations 
occurring in the numerical results. However, we can minimise these spurious 
oscillations by constructing a high-resolution scheme as discussed in Section 2.3, 
which also satisfies the TVD property. We use flux-limiter methods to adapt the 


Lax-Wendroff scheme for a system of conservation laws to a high-resolution scheme 
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by constructing a numerical flux-function of the form 


TVD _ FO SO FO 
FIV? = FEO + (PS — FES), (2.21) 
where ae is a second order numerical flux, BG is a first order numerical flux, 
2 

® = diag(®,) where ©, is a flux-limiter, which can be any of the flux-limiters listed 
in Table 2.1, and k is the jth component of the system. For the second order 
numerical flux, we use the Lax-Wendroff numerical flux (2.16) and re-write it by 
using 


A=XAX7}, 


where X is a matrix containing the right eigenvectors, e,, of A and A = diag(A,) is 


the diagonal matrix of eigenvalues, A, of A. Hence, we obtain 


LW 
Piya 


*y 


SS. 


fa FET) — 5 (KAX")2 (Fla — FY) 


S 


(Fay + EF?) ~ 5 


Ble ble 


(XAPX™ Pa (w7,., — w?’). 
For the first order approximation, we use the upwind numerical flux 


UP 
Pua 


1 nm n nm 
(Fea TF F?) _ 5 lAlins (why — w’') (2.22) 


le wl Re 


1 
(Pha + F?) - S(XIAIK YR (wh — Ww?) 


Thus, 


1 n n n 
Fa ans Le = 5 (X|A| = s|A|)X~") 541 (wii - w;’) 


and by using (2.6), we obtain the flux-limited Lax-Wendroff numerical flux 


1 n n n 
R= (X|A|LX~") 7) (w,,—wi), (2.23) 


(Fh + F*) ~ 5 


1 
2 
where 

L = diag (1 — ©(6,)(1 — s|al)) 
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A more convenient way of writing (2.23) is in a decomposed form, 


= 1 n 7: 1 ud . 
fg = 9 Pha + FY) — 5D lowlrell — 20%) — Mel enlieg, (2-24) 
k=1 
where 
(an) 7a . 
= (anh, », L=t= sen(vi)iy 15 ie. 
Here, p is the number of components of the system, k = 1,2,...,p represents the k*” 


component of the system and a, are the wave strengths, i.e. X~'Aw, associated with 
each component of the decomposition onto the eigenvectors e,. We can determine 


the values of a, by using 


3 


n 
= (nek) 542 where Aw, 1 = Wii — W; 


and for the shallow water equations, we obtain 


1 1 " 
Q19)",1 = = | Ah = —= (A(wh) — uAh : 

(analy = 5 [/ARF Tp Aluh) ~ wah . 
We require a source term approximation for the numerical scheme and adopt an 
approach discussed by Hubbard & Garcia-Navarro [18, 19], which is based on the 
approach of Roe [35] and Glaister [11]. The approach derives an upwind source term 
approximation that balances with the flux functions thus, making the numerical 
scheme satisfy the C-property. The source term approximation is constructed in a 
similar way to the numerical flux, 

R; = RU? +R, (2.25) 

2 2 

where 


TVD _ RFO SO _ PFO 
Riv _ Ria +® (RE Ri) , 
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For the first order approximation, we used the upwind flux (2.22) and adopt the 


approach of Bermudez & Vazquez [1] where 


Ble ble 


FO 
Riva 


(I+ |A|A™) R),.3 


(X (I= A™|A]) XR), 


pi 
2 


bie 


? 


which satisfies the C-property. Now, for the second order approximation, we used the 
Lax-Wendroff numerical flux. For the Lax-Wendroff flux, we deduced that a more 
accurate approximation of the source term was to use the second order accurate 
approximation (2.19), which satisfies the C-property. Thus, we use this second 
order approximation to balance the Lax-Wendroff numerical flux. Now, 


Ri“ = x=((1#s8A) R)j42 


NLR bw] rR 


= =(X(1¥sA7A’) X"'R) 


tt 


Nile 


Therefore ; 
Rit _ Rea =F, (XA~"|A| (s|A| — I) XR) 


its 


NI 


and by using (2.25) we obtain a flux-limited second order approximation of the 


source term 


= Ria + Rea, (2.26) 
2 2 

where 1 

SIE: at —1 —1 

na (X(I4A“|A|L) XR), 

and for the shallow water equations, 

m 0 

i+4 a g n n n n 

2 — 5 (hes + h?)(BE, — BP) 


Alternatively, we can re-write (2.26) in decomposed form, 


RE, = 5D [een( £ sam(As)(1— 8(04)(1— bal) 14 


1 
2 
k=1 


AT 


where the values of (3, are the components of X~'R and can be determined from 
P 
S> [6 pele = Ria ‘ 
k=1 


For the shallow water equations, we obtain 


(Bi2)ia = + 5 (va RAB) _ Where ABI = a 2B, 


and to satisfy the C-property, we must take a 5 (wr + wr). 


iw] 


It should be noted at this point that by applying flux-limiters to the source term 
approximation as well as the flux-function, spurious oscillations can still occur in 
the numerical results, especially if the source term is significant. For the shallow 
water equations, the source term becomes more significant as the variation of the 
riverbed becomes more pronounced. To minimise the chances of spurious oscillations 
occurring in the numerical results, we must choose the flux-limiter carefully. For 
example, the superbee flux-limiter reduces the CFL limit of the scheme more than 
the minmod flux-limiter and thus, the superbee flux-limiter has more of a chance of 
producing spurious oscillations than the minmod flux-limiter. Hudson [20] illustrated 
the importance of choosing the correct limiter and deduced that, in general, the 
minmod limiter was the most accurate limiter. In some cases, even with the minmod 
limiter, the scheme can still produce spurious oscillations but this is very rare and 


can usually be resolved by using a smaller Courant number. 


2.6.3  MacCormack Approach 


LeVeque & Yee [27] and Yee [45] adapted the MacCormack approach [28] to 


approximate systems of conservation laws with a source term present, 


1 8 
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where 


a 


wi) = w, — 3(Fy,, —F;) + sRjy1- 


The MacCormack scheme is a predictor-corrector scheme and reverts back to the 
classic Lax-Wendroff scheme for the constant coefficient case. The advantage of 
using (2.27) is that we do not need to approximate the eigenvalues and eigenvectors 
of the Jacobian matrix. However, the numerical scheme does not satisfy the TVD 
property, thus spurious oscillations may occur in the numerical solution. LeVeque 
& Yee [27] adapted (2.27) so that the numerical scheme satisfies the TVD property 
by using slope limiters 

wet = w?)4+D",—-Di (2.28) 

3 3 


where w.”) is the numerical approximation derived from (2.27), 


Dia = 2 > [lve ](1 — [vel (on — Qrenliga ’ 
and a, are the wave strengths associated with each component of the decomposition, 


i.e. X~tAw, and can be determined by using 
Pp 


n _ n _ 
Aw? = ) (nek) 542 where Aw; 41 = Wi41 — Wj. 
k=1 


The slope-limiter, (Qi). = Q* ,, can be any of the following (see Yee [45]) 
2 2 


i+ 
Qi = minmod(aj 1,074.1, v4.3) (2.29a) 
Qi = minmod(a7_1, oF 41) + minmod(a7, 1, 3) =O) 4 (2.29b) 
Qi = minmod(2a* 1, 20%, 1, 20%, a, sat y +a +3) (2.29c) 


where 


d min(|al ,|O|, el) if d = sgn(a) = sen(b) = sen(c), 


0 otherwise. 


minmod(a, b,c) = 
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For the shallow water equations, the source term approximation is 


0 
R” i ’ 
a |S (hts + h?)(Bha — BP) 


which makes the MacCormack scheme satisfy the C-property for the second order 
scheme (2.27). 


Unfortunately, the slope-limited version of the MacCormack scheme (2.28) does 
not satisfy the C-property, for consider the shallow water equations with the 
quiescent flow case, u = 0 and h = D— B. Here, the second term on the right 
hand side of (2.28) is 


dy + dy Ss 
D", = eee Am 4 gaa) gh oe O). 
ae Beas |. — oe eal 


Now, for the shallow water equations, 


1 
12 = 2 fan == (A(uh) _— uAh) 5 


1 
Vgh 
thus, for the quiescent flow case a, = sAh and so, Q; = Q2. Therefore, 


dy = dy = 5 / gh h(1 — s\/gh “Ah Q) 
and we obtain 


sVgh(1 — sV/gh)(sAh — Q) 


0 
i+$ 


n == 
Di4 > 


Hence, (2.28) does not satisfy the C-property unless Q = sAh, which reverts the 
scheme to the second order MacCormack approach. Unfortunately, the MacCormack 


scheme cannot satisfy the C-property by using slope-limiters but the scheme can if 
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we use flux-limiters instead. Hence, the flux-limited MacCormack approach is 


wrt! — w+D",—Dts (2.30) 
2 2 
where w.”) is the numerical approximation derived from (2.27), 
eat 
Di. = 3 [(a%m|Ax] — Pesgm(Ax)) (1 — [vel)( - ®(9x) en} 2 
k=1 


and © is now a flux-limiter, which can be any of the limiters listed in Table 2.1. As 
in the previous section, the values of 3; are the components of X~'R and can be 
determined from > 

> [eendia = Ry. 

k=1 


The source term approximation R’, must first be carefully chosen so that (2.27) 


dle 


satisfies the C-property. 


2.6.4 Numerical Results of the Different Lax-Wendroff 


Schemes 


To illustrate the accuracy of the Lax-Wendroff scheme and the flux-limited Lax- 
Wendroff scheme (2.23), we use Test Problem B, which is the wave propagation test 
problem, with w = 0.01. The Lax-Wendroff scheme (LxW) is used with both the 
pointwise (PW) source term approximation (2.18) and the second order accurate 
source term approximation (2.19), which satisfies the C-property (CP). The flux- 
limited (FxL) Lax-Wendroff scheme (2.23) is also used with the pointwise (PW) 
source term approximation (2.18) and the flux-limited source term approximation 
(2.26), which makes the scheme satisfy the C-property (CP). For both numerical 
schemes, we use v = 0.8, Ax = 0.01 and the minmod flux-limiter. A flux-limited 
second order version of Roe’s scheme, which is discussed in Section 2.7, on a fine 
mesh, Ax = 0.0001, is also used as a reference solution. The semi-implicit version 


of the schemes are not illustrated as the numerical results produced were practically 
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identical to the explicit schemes. From Figure 2.17 and Figure 2.18, we can see that 
both schemes with the pointwise source term approximation have produced poor 
numerical results due to spurious oscillations being present. The Lax-Wendroft 
scheme with the second order accurate source term approximation has produced 
considerably more accurate numerical results due to the numerical scheme satisfying 
the C-property but spurious oscillations are still present near the disturbance. The 
flux-limited scheme with the flux-limited source term approximation produced the 


most accurate numerical results due to the numerical scheme satisfying the C- 


property. 


2.6.5 Numerical Results of the Different MacCormack 
Approaches 


To illustrate the accuracy of the different MacCormack schemes, we use Test Problem 
C, which is the tidal wave propagation test problem. We use the second order (2.27), 
slope-limited (2.28) and flux-limited (2.30) versions of the MacCormack scheme. For 
all numerical schemes, we use v = 0.8, a spacial step-size of Ax = 6,480m and a 
final time of t = 10, 800s. For the slope-limited version of the MacCormack scheme, 
we use the limiter (2.29a) and for the flux-limited version of the MacCormack 
scheme, we use the minmod limiter. For this test problem, at t = 10,800s no 
movement should occur past x & 216,000m as the tidal wave has not propagated 
past this point. However, from Figure 2.19 and Figure 2.20 we can see that the 
slope-limited version of the MacCormack scheme has produced spurious waves past 
this point whereas the second order and flux-limited versions have not produced any 
movement past this point. This is due to the slope-limited version not satisfying the 
C-property whereas the other two versions satisfy this property. Hence, the slope- 
limited version is considerably less accurate than the second order and flux-limited 
versions of the MacCormack scheme. It should be noted that even though the second 
order version has produced very similar results to the flux-limited version, when a 
discontinuity is present the second order scheme will produce spurious oscillations 


due to the numerical scheme not satisfying the TVD property. 
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Figure 2.17: Numerical results of the different Lax-Wendroff schemes for Test 
Problem B with w = 0.01 at t = 0.7 (h+ B). 
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Figure 2.18: Numerical results of the different Lax-Wendroff schemes for Test 
Problem B with w = 0.01 at t = 0.7 (u). 
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Figure 2.19: Numerical results of the different MacCormack schemes for Test 
Problem C at t = 10,800 s (h+ B). 
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Figure 2.20: Numerical results of the different MacCormack schemes for Test 
Problem C at t = 10,800 s (u). 
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2.7 Roe’s Scheme 


Roe [33] derived an approach which approximates systems of conservation laws, 


Ow OF _ 


— 2.31 
a On (2.31) 


by using the numerical data, wi’, to construct a piecewise constant function 
n\ __ n n n 
wi(z,t")=wy, rE Eure aee 


Now, for any two adjacent states of the piecewise data, wy and wp, (2.31) can be 


solved by determining the exact solution of a linearised Riemann problem, 


Ow ~ Ow 

—+A — =(0 

gp ee) 

which is related to (2.31) where A(w ,, wr) © gE is the linearised Jacobian matrix 
and ~ is called the Roe average. The Roe averaged Jacobian matrix, A(w L, Wr), 


must be chosen such that the following properties are satisfied 


e A(w,, wp) must be diagonalisable with real eigenvalues (hyperbolicity); 


e A(wr,wr) — A(w) as wL,wR — w (consistency); 


e AF = A(w_z, wr)Aw (conservation). 


The eigenvalues and eigenvectors of A are \ and é respectively which are determined 
from the decomposition 
P 
k=1 
where Aw = wr - wy, p is the number of components in the system and da, 


represents the wave strengths associated with each component of the decomposition, 
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i.e. X~!Aw, and are determined from 


D 
k=1 


Once the eigenvalues, eigenvectors and wave strengths associated with the linearised 
Riemann problem have been obtained, we can use any of the numerical schemes 
discussed in the previous section with the Roe averaged values. For example, the 


flux-limited version of the Lax-Wendroff scheme can be used, 


. 1 i - 
Fi =5at+F-5>> |e al( (1 @(0,)(1 — rel) @e] (2.32) 
k=1 v2 
where 
~ (Qk) r44 ; 
Ve = SAn, On = Cee [=1- sgn (Ve) i421, 
1? 


and ®; can be any of the flux-limiters listed in Table 2.1. Since this version of 
the flux-limited Lax-Wendroff scheme uses the Roe averaged values, we call this a 
flux-limited version of Roe’s Scheme. Glaister [11, 12] determined the following Roe 


averaged Jacobian matrix for the shallow water equations 


with eigenvalues 


peal 
Cy 
Or 
iY) 
=| 
Qa. 
pez 
i) 

| 
oy 
+ 
Or 


and corresponding eigenvectors 


&=|_ ._ and @» ae 
u—=C UtC 
and wave strengths 
1 i 
Q1.2 = gah 2 (A(uh) = uAh) 
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where the Roe averages are 


Vhrur + Vhcur 
VhatVht 


The Roe averaged values can also be used for the source term approximations. For 


m= 


and h= 5 (hehe). 


example, the source term approximation that was derived for the flux-limited Lax- 
Wendroff scheme (2.26) can be used, 


f= .> | Bex + sen(Ax)(1 — ® (1 — [rel))) (2.33) 


1 . 1 
k= 2 


where the values of Br are the components of X~!R and are determined from 
1G; s 
i > Gre, = R. 
Glaister [11, 12] also determined that for the shallow water equations, 
~ CAB cAB 


A= ; and $y = — - 


2.7.1 Numerical Results of the Flux-Limited Version of Roe’s 


Scheme 


To compare the accuracy of the Roe averaged Jacobian matrix with the averaged 
Jacobian matrix (2.17), we use the flux-limited Lax-Wendroff numerical flux (2.23) 
with the source term approximation (2.33), which satisfies the C-property. Test 
Problem C is used with a final time of t = 21,600s and both schemes are run 
using vy = 0.9, a spatial step-size Ax = 6,480m and the minmod flux-limiter. 
Figure 2.21 illustrates a comparison of the numerical results of the surface 
elevation for both schemes. Here, we can see that both approximations have 
produced very similar numerical results. The numerical results are smooth and 


are free from spurious oscillations. Both schemes have produced similar numerical 


57 


70 


Surface elevation 


a 


60 


10) 108000 216000 324000 432000 540000 648000 


x 


—& Averaged —— Roe Averaged 


Figure 2.21: Numerical results of the flux-limited scheme with the different Jacobian 
approximations for Test Problem C at t = 21,600s (h + B). 


results due to the approximations of the Jacobian matrices being similar. For the 
shallow water equations, the only difference between the two approximations of the 


Jacobian is the approximation of u, 


e = Vhrur + Vhypuy 
VhrtVht 


as the approximation of h is the same. This results in both approximations of 
the Jacobian matrices being similar for most test problems. However, the two 
approximations of the Jacobian matrices will become increasingly different as the 


system grows due to the complexity of the Jacobian. 


58 


2.8 Summary 


Throughout this chapter, we have discussed a variety of numerical scheme that 
can be used to approximate systems of conservation laws with source terms. We 
have shown that the LxF and NT schemes are accurate when a source term is 
not present but inaccurate when a source term is present due to the numerical 
schemes not satisfying the C-property. We have derived a flux-limited version of 
the Lax-Wendroff scheme, which satisfies the C-property, and is very accurate even 
when the source term is significant. A flux-limited version of the MacCormack 
scheme was also derived that is also accurate for all of the test problems discussed 
in this chapter. Two approximations of the Jacobian matrix, a basic averaged 
approximation and a Roe averaged approximation, were derived for the flux-limited 
Lax-Wendroff scheme. For the shallow water equations, both approximations of 
the Jacobian matrix produced very similar results. However, for larger systems the 
Jacobian matrix becomes more complex, which may result in the two approximations 
of the Jacobian matrix producing different results. In the next chapter, we extend 
the schemes discussed in this chapter to approximate the equations governing 


sediment transport in one dimension. 
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Chapter 3 


Numerical Formulations for 
Approximating the Equations 
Governing Sediment Transport in 


One Dimension 


In the previous chapter we discussed a variety of numerical schemes that can be used 
to approximate systems of conservation laws in one dimension. We used the shallow 
water equations to illustrate the derivation and accuracy of the different numerical 
schemes. In this chapter we adapt the numerical schemes discussed in the previous 
chapter to numerically approximate the equations governing sediment transport in 
1D. The classic Lax-Wendroff scheme and the flux-limited version of Roe’s scheme 
are adapted to approximate the equations. The classic Lax-Wendroff scheme is 
widely used in industry, but not surprisingly the numerical results obtained suffer 
from spurious oscillations resulting in the numerical scheme becoming unstable for 
long time periods. Different measures have been applied to the classic Lax-Wendroff 
scheme to try to eliminate the spurious oscillations such as a smaller Courant number 
and the scheme has been adapted to satisfy the TVD property. Unfortunately, the 


spurious oscillations could not be eliminated and overpowered the numerical results 
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for long computational run times, see Damgaard [4] and Damgaard & Chesher [5]. 
These difficulties are mainly due to the source term approximation, i.e. the scheme 
does not satisfy the C-property. The flux-limited version of Roe’s scheme was derived 
to overcome these difficulties so that we can obtain accurate numerical results with 
no spurious oscillations present, even for long time periods. Thus, it is useful to 
illustrate the numerical results obtained from both schemes to highlight the improved 
accuracy of the flux-limited version of Roe’s scheme compared to the classic Lax- 
Wendroff scheme. 


Before we adapt the numerical schemes, we look at five different formulations 


that can be used with the numerical schemes to approximate the equations. 


3.1 Different Formulations 


There are numerous approaches that can be used to obtain an approximation of the 
equations governing sediment transport. In this thesis, we consider following two 


approaches: 


1. the steady approach, where the water flow is assumed to be steady and the 
changes in the bed update have a negligible effect on the water flow, i.e. the 
wave speed of the bed-updating equation is considerably smaller in magnitude 
than the wave speeds of the water flow. By making these assumptions, the 
system is decoupled into a water flow approximation, which is iterated to an 


equilibrium state, followed by a bed update. 


2. the unsteady approach, where no assumptions are made and the water flow 
and riverbed are calculated simultaneously. With this approach, the water 
flow can be either steady or unsteady and the changes in the bed update are 
considered to be significant, i.e. the wave speed of the bed-updating equation is 
of a similar magnitude to the wave speeds of the water flow. For this approach, 


the system is discretised simultaneously. 
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A considerable amount of research has been carried out on the steady approach, 
which was pioneered by Cunge et al. [3] who discussed both approaches at length. 
Cunge et al. stated that for most physical cases, the bed moves at a considerably 
slower speed than the water flow. For example, for a 100 km channel, the water 
flow usually takes one or two days to propagate over the entire channel whereas the 
bed usually takes many years. Thus, we expect the wave speed of the bed-updating 
equation to be considerably smaller in magnitude than the wave speeds associated 
with the water flow. This enables the steady approach to be derived, where the 
water flow and the bed update is decoupled and discretised separately making the 
equations considerably easier to numerically approximate. The steady approach of 
Cunge et al. also assumes the water flow is in an equilibrium state which allows 


equations (1.1) and (1.2) to be re-written as 


Qe =0 => Q(a, t) = Qe V(z, t), (3.1) 


which results in a constant discharge throughout the domain, and 


pe + g(h+ a) =0, (3.2) 
where the equation has been re-written in nonconservative variable form. Thus, 
a numerical approximation of (3.2) and (1.3) is now only required. In this thesis, 
we slightly adapt the approach by iterating the water flow to an equilibrium state, 
which has the effect of imposing a constant discharge, and then updating the bed. 
Unfortunately, the steady approach is very limited as the water flow is assumed to 
be in an equilibrium state, but rarely is. There are also cases, where the wave speed 
of the bed-updating equation is of a similar magnitude to the wave speeds of the 
water flow, see Perdreau and Cunge [31], but this case is rare. However, the steady 
approach cannot be used for these cases as the bed is assumed to have a negligible 


effect on the water flow. 
The unsteady approach approximates the whole system simultaneously on the 


same time scale. This approach can be used with steady and unsteady water flow 


and is considerably more robust than the steady approach. Little research has been 
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carried out on the unsteady approach as until recently, a full numerical solution of 
the approach was considered to be too complex and expensive to obtain, see Cunge 
et al. [3]. This is mainly due to the complexity of the sediment transport flux, 
which can range from a basic analytical function, as discussed in Section 1.4, to a 
system of equations such as the equations derived by Einstein [8]. In some cases, 
the sediment transport flux cannot be written analytically and is determined by 
using a “black box” approach where the flux is deduced from empirical data. The 
steady approach can easily incorporate all of these sediment transport fluxes, but 


the unsteady approach is considerably harder. 


In this thesis, to determine the validity of the steady and unsteady approach, we 
consider five different formulations that can be used to approximate the equations 
governing sediment transport. Each formulation will be derived for the sediment 


transport flux (1.10) for positive values of u, 
q(u) = Au” where u> 0. (3:3) 


Four of the formulations will be based on the unsteady approach and one will be 


based on the steady approach. 


3.1.1 Formulation A 
The equations (1.1), (1.2) and (1.3) can be used as written in three different ways: 


1. Formulation A-CV: The following method is based on the steady approach 


and consists of 


e a “fixed-bottom step”, where the shallow water equations, 


Beh 


are iterated to an equilibrium state whilst keeping the bed fixed. 


0 
—ghB, 


pI 


uh 7 
hu? + igh? ; = 
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e a “changing bottom step”, where the bed is updated whist keeping all 


other variables fixed. 


Numerically, an equilibrium state has been obtained when 
wet? —w?| <tol Vi, 


where tol is the desired tolerance level. Thus the numerical approximation 
of the water flow must satisfy the tolerance level before the bed is updated 
and the process is repeated. The overall time step of this formulation is the 


morphological time step of the bed-updating equation. 


. Formulation A-NC: Here, the shallow water equations and the bed-updating 
equation are numerically approximated sequentially using the same time step. 
This formulation is similar to Formulation A-CV as the water flow is still 
calculated separately from the bed update. However, the water flow is no 


longer iterated to an equilibrium or steady state after each bed update. 


. Formulation A-SF: All three equations are re-written in system form, 


h uh 0 
uh | +] hu?+$gh? | =| —ghB, (3.4) 
B |, &q . 0 


and the whole system is numerically approximated simultaneously. This 


formulation is based on the unsteady approach. 


All three formulations have a source term present, which must be treated carefully 


to avoid difficulties obtaining an accurate numerical approximation. The Jacobian 


matrix will be required for each variation of Formulation A. For 


1. Formulations A-NC and A-CV, the shallow water equations’ Jacobian matrix 


is 


where c = /gh, which has eigenvalues 
Ay=u-—c and A»=utc, 
with corresponding eigenvectors 


| 


t= 6 


1 
U+tC 


e, = and e)= 


Both Formulations A-CV and A-NC require an approximation of the wave 


speed of the bed-updating equation, 
Og 


which can be difficult to obtain as the sediment transport flux is not a direct 


function of B, see Section 3.4. 


. Formulation A-SF, if the sediment transport flux (3.3) is used then the 


Jacobian matrix is 


0 1 O 
A=| e=e? 2u 0 |, 
—ud d 0 


where c = /gh and d= £ Amu”! for u > 0. Notice that this Jacobian matrix 
is singular, which we might expect to create difficulties when implementing a 
numerical scheme for this formulation. The eigenvalues of the Jacobian matrix 
are 


Ay=u-—c, Ax=O0 and A3=ut+c, 


with corresponding eigenvectors 


1 0 1 
a=] Y-C], eae =| 0 and eg3= | ute 
—cd 1 cd 
U—C utc 


For this formulation, notice that the eigenvalue associated with the bed- 
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updating equation is zero due to the Jacobian being singular. This implies 
that movement in the bed is created only by the water flow and thus, implies 


that the water flow cannot be decoupled from the bed update. 


3.1.2 Formulation B 


Another approach that can be used is to re-write the equation of conservation of 


momentum (1.2) as 
du  O[fu?+g(h+ B)| 


oo a = 0, (3.5) 


by using (1.1) and then combine (1.1), (3.5) and (1.3) into system form to obtain 


Formulation B, 


h uh 
u | + | gu?+g(h+B) | =0. (3.6) 
B t 6q x 


Notice that this formulation does not have a source term present, hence, Formulation 
B will be easier to approximate numerically. However, Formulation B may not be 
in conservative variable form, which could results in shocks propagating at incorrect 
speeds. If the sediment transport flux (3.3) is used then the Jacobian matrix of 


Formulation B is 
u h O 


A=|/g ug], 
0d 0 


where d = A€mu'~' for u > 0. Notice that the Jacobian matrix of Formulation B 
is not singular. The eigenvalues, A, of the Jacobian matrix cannot be easily written 


analytically since they are the roots of the polynomial 


P(A, w) = 3 — 2ud? + [u? —g(h+ d)| A+ gud = 0. 
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However, we can prove that the roots of P(\, w) are always real by using formulae 


for the roots of a cubic, see Spiegel & Liu [37]. For a cubic equation, 
P(x) = 2? + a,x" + agx + az = 0, 


if we let 


I 1 
= g (a2 —a?) and R= 5g (90182 — 27a3 — 2a‘), 


then the discriminant is D = Q? + R? and if 


1. D> 0 then one root is real and two are complex; 
2. D=0 then all roots are real and two are equal; 


3. D <0 then all roots are real and unequal. 


If D <0 then the roots of P(x) are 


1 1 
1 = 2\/—Q cos(=8) — =a, (3.7a) 
3 3 
1 1 
Lo = 2\/-Q cos(3 (0 + 2m)) — 3 (3.7b) 
and : ; 
£3 = 2,/-Q cos( (0 + 47)) — gu (3.7c) 
where cos @ = 2 


o/ —Q 
Now, by using the above approach, for Formulation B 
a, =—2u, ag=u?—g(h+d) and a3=gud 


which implies that 


U 


5q(99(2h — d) — 2u°). 


Q=-J(w+39(h+d)) and R= 


67 


Hence, 
De FalBauh? — u?(Ahu? + gd(d + 20h)) — 4g?(h3 + d? + 3hd(h + d))| 
and for all three roots to be real and unequal, 


8gu7h? < u?(4hu? + gd(d + 20h)) + 4g°(h® + d® + 3hd(h + d)), 


which is always satisfied since h(x,t) > 0 and u(x,t) > 0 > d> 0. Hence, the roots 
of P(A, w) are always real and unequal. The signs of the roots of the Jacobian can 


be obtained from the determinant of A 


and since, the determinant is negative this implies that either all three eigenvalues 
are negative or two are positive and one is negative. Also, as the parameter A — 0, 


the eigenvalues of the Jacobian tend to 
Ayau-c, Ag>0 and A3—> ute, 


which are identical to the eigenvalues of Formulation A-SF. Thus, we can determine 


that one of the roots is negative and two of the roots are positive, i.e. 
M1 <0< r2 < A3. 


De Vries [7| obtained an approximation of the eigenvalue associated with the bed- 


updating equation by assuming 
Ay=u-—c and A,s=utc, 


and substituting into (3.8) to obtain 
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De Vries deduced from this approximation that A, < 0 and A» > 0 for subcritical 
flow, wu < c, whereas A; > 0 and A2 < 0 for supercritical flow, u > c. Now that we 
have proved that the roots of P(A) are always real, the eigenvectors, 
1 
1 
(u— Ax)? — gh 
gh 


can be obtained where the A; are given by (3.7). 


3.1.3. Formulation C 


In order to obtain a formulation that is written in conservative variable form and 
whose Jacobian matrix is not singular, we use the product rule, (hB), = hB,+Bhz, 
to re-write the equation of the conservation of momentum (1.2) as 

O(uh)  O[hu? + dgh? + ghB] 


en a a =i: (3.9) 


Now, combining equations (1.1), (3.9) and (1.3) into system form, we obtain 


h uh 0 
uh | + | hu?+3gh?+ghB | =| gBhz |. (3.10) 
B I &q . 0 


Unfortunately, Formulation C has a source term present. If the sediment transport 


flux (3.3) is used then the Jacobian matrix of Formulation C is 


0 1 O 
A= | g(h+B)—wu? 2u gh], 
—ud d O 
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where d = £ Amu”! for u > 0. The eigenvalues, A, of the Jacobian matrix again 


cannot be easily written analytically since they are the roots of the polynomial 
P(A, w) = »° — 2ud? + [u? — g(h + B+ hd)] A+ ghud = 0. 


However, we can prove that the roots of P(\,w) are always real by using the 


approach discussed in Section 3.1.2. For Formulation C, 
@,)=—-2u, ag=u?—g(ih+B+hd) and ag =ghud 
which implies that 


1 
Q=-F(w + 39(h+B+hd)) and R= — 5 (2u" — 9g(2h + 2B — hd)). 


Hence, 


D= Fa lbura(h +B)? — wW[4u2(B + h) + ghd(20(h + B) + ha)] 


—497(h? + Be + hid? + 3A(d+1)(h+ B)(hd+ B))| 


and for all three roots to be real and unequal, 


8u29(h +B)? < u?[4u?(B +h) + ghd(20(h + B) + hd) 
+497(h? + BS + h'd? + 3h(d+1)(h+ B)(hd+ B)), 


which is satisfied if h(xz,t) + B(a,t) > 0 since h(z,t) > 0, u(z,t) > 0 > d>0 and 
—oo < B(x,t) < oo. Hence, the roots of P(A, w) are always real and unequal if 
h(a,t)+ B(z,t) > 0. As with Formulation B, it can also be shown that P(X, w) has 


one negative root and two positive roots, i.e. 


My <0 < Aq < Az. 


70 


Now that we have proved that the roots of P(\) are always real, we need to determine 


the eigenvectors. The eigenvectors of the Jacobian matrix are 


1 
rx 

u? — g(h + B) + (Ap — 2u)Ax 
gh 


er = 


? 


where again A, are given by (3.7). 


Throughout this section, we have discussed five different formulations that can 
be used to numerically approximate the system (1.1), (1.2) and (1.3). In the next 
section, we discuss how these different formulations can be numerically approximated 
using the flux-limited version of Roe’s scheme, which was discussed in the previous 


chapter. 


3.2 Adaptation of the Classic Lax-Wendroff 


Scheme 


The classic Lax-Wendroff scheme (2.16) can be used to approximate the equations 
governing sediment transport. Unfortunately, the scheme suffers from spurious 
oscillations, which can result in an inaccurate numerical solution. For the source 
term approximation, we adopt the pointwise approach, which is used in industry. 
The classic Lax-Wendroff scheme does not satisfy the C-property when the pointwise 


source term approximation is used. 


The classic Lax-Wendroff scheme can easily be adapted to approximate the 


different formulations discussed in the previous section. 


rae 


3.2.1 Formulation A-CV 


For Formulation A-CV, the system is decoupled into a water flow approximation 
followed by a bed update. The classic Lax-Wendroff scheme can easily be adapted to 
approximate this formulation. In the previous chapter, the shallow water equations 
were approximated using the classic Lax-Wendroff scheme, see Section 2.6.1. 
However, for the bed-updating equation we need to adapt the scheme to numerically 


approximate the scalar conservation law. For the bed-updating equation, we obtain 
Bit! = BP —€s (qt, —at 3), (3.11) 
2 2 
where the numerical flux-function is 
a4. — 9 (Cam + qi’) = att (Ca —- qi) 


and the wave speed approximation is 


dq |" 

- OB ith 
Unfortunately, (3.11) requires an accurate approximation of the wave speed, which 
can be very difficult to obtain since the sediment transport flux is not a direct 
function of B. A variety of approaches that can be used to obtain an approximation 
of the wave speed are discussed in Section 3.4 where the advantages and 


disadvantages of each approach are discussed as well. 


3.2.2 Formulation A-NC 


Formulation A-NC is practically identical to Formulation A-CV, with the exception 
that the shallow water equations are no longer iterated to an equilibrium or steady 


state. 
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3.2.3. Formulation A-SF 


The classic Lax-Wendroff scheme (2.16) can be used to approximate Formulation 
A-SF as written. The Jacobian matrix is given in Section 3.1 and the source term 


is approximated using a pointwise approach, 


0 
R; a —$ (hh + he 1)( Bry ~~ B71) 
0 


3.2.4 Formulation B 


The classic Lax-Wendroff scheme (2.16) can be used to approximate Formulation 
B. The Jacobian matrix is given in Section 3.1 and no source term is present for 
Formulation B, thus 

RR = 0. 


4 


3.2.5 Formulation C 


The classic Lax-Wendroff scheme (2.16) can be used to approximate Formulation 
C as written. The Jacobian matrix is given in Section 3.1 and the source term is 


approximated using a pointwise approach, 


0 
Rj = (Bia + BE) (hE — i) 
0 
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3.3 Adaptation of the Flux-Limited Version of 


Roe’s Scheme 


In the previous chapter, we discussed a variety of numerical schemes that can be 
used to approximate systems of conservation laws with source term. We used the 
shallow water equations to illustrate our techniques and to determine which scheme 
produced the most accurate numerical results. We determined that out of all of the 
numerical schemes discussed, the flux-limited version of Roe’s scheme (2.32) with the 
source term approximation (2.33) was the most accurate numerical scheme. This was 
mainly due to the scheme satisfying the C-property. Thus, in this section we adapt 
the flux-limited version of Roe’s scheme to approximate the different formulations 


discussed in the previous section. 


3.3.1 Formulation A-CV 


For Formulation A-CV, the system is decoupled into a water flow approximation, 
the shallow water equations, and a bed update. The flux-limited version of Roe’s 
scheme can easily adapted to approximate this formulation. In the previous chapter, 
this scheme was discussed based on the shallow water equations. However, for bed- 
updating equation, we need to adapt the flux-limited version of Roe’s scheme to 
numerically approximate the scalar conservation law. For the bed-updating equation, 
we obtain 


Bp = Bt — 6 (ah, a 


) (3.12) 


i 
2 


where the numerical flux-function is 


n 


1 
qiy = 3 (ata?) -—= 


D 5 [Nea (1 © (is) (2— fey) (Bhs — BP), 
n n n Bry = BY 
a SArt,  OFa BY, — BY 


74 


F=imsen (viiy), y= esa] 

i+} 
and ® can be any of the flux-limiters listed in Table 2.1. In order to be able to 
use (3.12), we need to be able to obtain an accurate approximation of the wave 
speed, A, which can be very difficult to obtain since the sediment transport flux is 
not a direct function of B. A variety of approaches that can be used to obtain an 
approximation of the wave speed are discussed in Section 3.4 where the advantages 


and disadvantages of each approach are discussed as well. 


3.3.2. Formulation A-NC 


Formulation A-NC is practically identical to Formulation A-CV, with the exception 
that the shallow water equations are no longer iterated to an equilibrium or steady 


state. 


3.3.3. Formulation A-SF 


Formulation A-SF can easily be used with the flux-limited version of Roe’s scheme. 


For the sediment transport flux (3.3), we obtain the Roe averaged Jacobian matrix 


0 1 0 
A=|@-@ 2% 0 |, 
tid d 0 


where ¢ = \/ gh. The eigenvalues of the Roe averaged Jacobian matrix are 


M=t—-G M=0 and \3=t4+e 
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with corresponding eigenvectors 


1 0 1 
e, = we 5 eo = 0 and e3 = ue 
—éd 1 cd 
u—¢ Ute 

and the wave strengths are 

= LAR . (A(uh) — uAh) 

“9 gg 

i — (a? + @)Ah 
ip = AB — A oe ) 
(i — ¢)(u+ 6) 

and i i 

a3 = gah + 2G (A(uh) — uAh). 
The Roe averaged values are 

Vvh +WVh oe | 

_= ee hee 
VhaetwVht 2 
and 
PS EA(Au™) 
~ A(uh) — aAh 
If m is an integer, then d can be re-written as 
a A h m-1 
A= E(Vhrt vhz) (uz) ay, 
For the source term approximation, 
2 CAB x die? AB P cAB 
_ SS d == ‘ 
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3.3.4 Formulation B 


Formulation B can also be used with the flux-limited version of Roe’s scheme. For 
the sediment transport flux (3.3), we can obtain the following Roe averaged Jacobian 


matrix 


> 


0 


g 
0 


A= 


one f& 
Quy 2 


The eigenvalues, \, of the Roe averaged Jacobian are obtained by finding the roots 


of the polynomial 
P(X) = 33 — 232 + ji? —g(h+d)| \+ gid =0. 


The roots of P(A) are determined by using the approach which was discussed in 
Section 3.1.2. Once the Roe averaged eigenvalues have been obtained, they are used 


to determine the values of the corresponding eigenvectors, 


and the wave strengths 


k= Na) (& — ds) + gh) AR + (2% — Ag — Ax)/hAu+ ghAB 
k —= =. ~ ~ ~ 
(An — Aa) (Ax — As) 


where a #4 k £b. The Roe averaged values are 


EA(Au™) 
Au 


1 ~ 1 7 
i= 5 (ur + uz), h= 5(hr + hz) and d= 
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If m is an integer, then d can be re-written as 


Since Formulation B does not have a source term present, 


3.3.5 Formulation C 


We can also adapt the flux-limited version of Roe’s scheme to approximate 


Formulation C. For the sediment transport flux (3.3), the Jacobian matrix is 


0 1 0 
A=| g(h+B)—-w 2% gh 
—iid d 0 


The eigenvalues, , of the Roe averaged Jacobian matrix are obtained by finding 


the roots of the polynomial 
P(A) = X¥ = 2? + |i? — g(h + B+ hd)] A+ ghiid = 0. 


The roots of P(X) are determined by using the approach which was discussed in 
Section 3.1.2. Once the Roe averaged eigenvalues have been obtained, they are used 


to determine the values of the corresponding eigenvectors, 


1 
Ne 

a? — g(h+ B) + (Ap — 28) rg 
gh 


e, = 
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and wave strengths 


(gdp + g(h + B) — @)Ah + (2% — Aq — Ax) A(uh) + ghAB 
Oe — Aa)(Ae — Av) 


Ak = 


where a#k #b. The Roe averages are 


V hpur ar AL hypuz 
Vhr+ Vht 


- <4 si 
> h= 5 lhe t hr), B= 5(Br+ Br) 


u= 


and 
EA(Au™) 


= A(uh) — GAR’ 


If m is an integer, then d can be re-written as 


1 


AE(Vha + Vhz) © 


d= 
Vv hyhr sy hrhy =, 


(ur)*(uz)™ Oo. 


For the source term approximation, we use 


Be = = —-- %) where a#Ak+£b. 


3.4 Approximating the Wave Speed of the Bed- 
Updating Equation 


Most of the numerical schemes that can be used with Formulations A-NC and A-CV 


require a numerical approximation of the wave speed of the bed-updating equation, 


whose wave speed is A = & oe. Unfortunately, the sediment transport flux, q, is 
not a direct function of B, which can create difficulties in obtaining an accurate 


approximation of the wave speed. 
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q(u) 


Figure 3.1: The occurrence of a negative wave speed. 


One approach that can be used to approximate the wave speed is to use a finite 
difference approach, 


Gi+1 — Gi 


Unfortunately, the finite difference approach can only be used when Bj, — Bi 4 0 
and can also produce an inaccurate approximation of the wave speed when the 
gradient of the riverbed changes sign. If the gradient of the riverbed and the velocity 
do not change sign at the same time, the sign of the wave speed may change. For 
example, for a small positive velocity, the wave speed of the bed-updating equation 
is always positive. However, the finite difference approach may produce a negative 
wave speed, as illustrated in Figure 3.1, where it can be seen that Bj., — B; > 0 but 


ditt — G <0, which implies that for the illustrated data 


\dit1 _ qi| 
re ee oe 
b= Sop By 


1 
2 


Hence, the finite difference approach has produced a negative wave speed but the 


analytic wave speed is positive. 


Another approach that can be used is an analytical approach as discussed by 
Chesher et al. [2], but can only be used for problems where the height of the water 
from the bottom of the channel is considerably larger than the height of the riverbed. 
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The approach assumes the river is of a constant height, h +B» D, and Q = uh so 


that u can be re-written in terms of B, 
u=Qh'sQ(D-B)t. 


As an example, consider the sediment transport flux (3.3), which can be re-written 
in terms of B as 
q(u) = Au™ = AQ™(D — B)-™. 


Thus, the sediment transport flux can now be differentiated with respect to B, 


oq 
a © AmQ™(D-— B)™*. 
ot x AmQ"(D ~ B) 
Hence, we obtain an analytical approximation of the wave speed for the bed-updating 
equation, 
Am 
FY sau (3.14) 
This approximation of the wave speed for the bed-updating equation is very limited 
as it can only be used for subcritical flow, u < gh. For supercritical flow, the 
bed should propagate upstream but the analytical approximation of the wave speed 


results in the bed propagating downstream. 


De Vries [7] also derived an analytical approximation of the wave speed associated 
with the bed-updating equation, 
ud O 
AQ = g ; where ane [54], 


ce —u Ou 


and was discussed in Section 3.1.2. Unfortunately, the wave speed approximation is 
only valid for sediment transport fluxes that are know functions of wu only, i.e. q(w), 
and the approximation cannot be used when the flow is near critical. However, the 
approximation can be used for subcritical and supercritical flow. For the sediment 
transport flux (3.3), d= Amu™! thus, 

_ A€mgu™ 
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Wave speed approximation 


Oo 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
Fr 


—&- Chesher et al —— De Vries 


Figure 3.2: Comparison of the analytical wave speed approximations (3.14) and 
(3.15) with D = 10 and A= 0.001 for F, = 0 to 2. 


From Figure 3.2, we can see that the analytical wave speed approximations of De 
Vries and Chesher et al. produce similar values when the Froude number is small, 
ie. F. < 0.4. In this thesis, only small Froude numbers are considered thus, the 


wave speed approximation derived by Chesher et al. (3.14) is used. 


3.5 Test Problems 


To determine which formulation with the flux-limited version of Roe’s scheme 
produces the most accurate numerical results, we use the following two test problems. 
Both test problems consist of a channel of length 1000m with the following dummy 


initial conditions 


Q 


h*(@,0) =10— BY(e,0) and u*(z,0) = Ep 
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where @ is a constant. The bathymetry for 


1. Channel Test Problem A is 


Bein? (76t= 300) 
200 


0 otherwise 


if 300 < x < 500 
B(z,0) = 


r) 


where B is the maximum height of the bed. 


2. Channel Test Problem B is 


B(c,0) By we 300 
20) = : 
Br if x > 300 


To obtain the initial conditions, the water flow is iterated with the dummy initial 


conditions to an equilibrium state whilst keeping the riverbed fixed, i.e. 
[wit —w?'| < tol, 


where tol is the desired tolerance level. This ensures that the initial conditions 
of the water flow and bed are consistent and greatly reduces the possibility of an 
impulsive start occurring for larger values of A. The initial conditions for Channel 
Test Problem A with B = lm and Q = 10 are illustrated in Figure 3.3 and Figure 3.4 
and for Channel Test Problem B with By, = 1m, Be = 0m and Q = 10 are illustrated 
in Figure 3.5 and Figure 3.6. 


For both test problems, the sediment transport flux discussed by Grass (1.10) is 


used with m = 3, 


which is valid for all values of u. To ensure the error of the numerical schemes do 


not grow, the variables are non-dimensionalised, 


B ae 
See 


Velocity 


Po | 


100 200 300 400 500 600 700 800 900 
x 


—E- B(x,0) —<— h*(x,0) + B*(x,0) —+— h(x,0) + B(x,0) 


Figure 3.3: Initial conditions for Channel Test Problem A (h & B). 


1000 


0 100 200 300 400 500 600 700 800 900 


—<— u*(x,0) —+— u(x,0) 


Figure 3.4: Initial conditions for Channel Test Problem A (uw). 
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1000 


10 * >< OOO a. oo oo oS 


9.998 4 


ys ee 


Velocity 


400 500 600 700 800 900 1000 
x 


—E- B(x,0) —<— h*(x,0) + B*(x,0) —+—h(x,0) + B(x,0) 


Figure 3.5: Initial conditions for Channel Test Problem B (h & B). 


Oo 100 200 300 400 500 600 700 800 900 1000 


—><— u*(x,0) —+— u(x,0) 


Figure 3.6: Initial conditions for Channel Test Problem B (uw). 
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, Ue ,» _ AL q tolt = te 
wes =F_ and tol'=—, 


where L = 1000 and T' = = are the non-dimensional coefficients. Here, L is 


taken as the length of the channel. 


It should be noted that even though in Chapter 1 it was stated that the differential 
form of equation for conservation of momentum (1.2) becomes invalid when a 
discontinuity appears in the riverbed, this is only analytically and the equation holds 
numerically as long as the discontinuity in the riverbed is small. Thus, Channel Test 


Problem B can be used as long as the discontinuity in the riverbed is small. 


3.5.1 Approximate Solution for Channel Test Problem A 


We can determine an approximate solution of Channel Test Problem A by assuming 
that the total height of the river is constant and the discharge is constant throughout 


the whole domain, i.e. 
hig.) 2D = Bet) and. ule,tih(z,t) = Qe -Vle,t), 


where Q, is a constant. These assumptions are only valid when the riverbed is 
interacting slowly with the water flow, A < 0.01, and when the water flow is moving 
slowly, Q. < 10. However, the assumptions enable the velocity to be re-written in 


terms of the riverbed, 


An approximate solution of B(x,t) is now required. The bed-updating equation 


(1.3) can be re-written in quasi-linear form as 


B, + AB, = 0, 
where the wave speed is 
Oq 
a Ea 


Since the sediment transport flux is not a direct function of B, the wave speed of 
the bed-updating equation can be difficult to obtain. However in Section 3.4, an 
analytical approximation was derived which assumes the total height of the river is 
constant and uses u * Q(D — B)~! to re-write the sediment transport flux (3.3) in 


terms of the riverbed, 
q(u) = Au™ = AQ™(D— B)™. 


The sediment transport flux can now be differentiated with respect to B and we 


obtain the analytical approximation of the wave speed 
As AEmQ™(D-By™), 


By assuming the discharge is constant, the characteristics are given by 


& = AgmQMD ~ B(x, 0)" 


and with initial bathymetry 


A — 300 
Bein? (™2— 30) it 300 < x < 500 
B(x,0) = 200 
0 otherwise 
thus 
—(m+1) 
: — 300) 
d _ +2 (Xo s <= < 
= = AémQ” (v Bsin ee if 300 < x < 500, 
D-(™+1) otherwise. 


Now, by integrating we obtain 


re . —(m+1) 
(D — Beir? (22553) ) if 300 < xq < 500, 


200 


t= to Asn, t 
Dower) otherwise, 


(3.16) 
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which unfortunately cannot be re-written in terms of x7 9. Hence, the approximate 


solution of B is 


Bsin? 
200 


0 otherwise, 


: mo) if 300 < ao < 500, 
Bet = 


where the value of x is determined by substituting values of xo and t into (3.16). 
Unfortunately, this approximate solution is only valid until the characteristics cross, 
which first occurs in the region 300 < xg < 500. We can determine the time the 


characteristics first cross by using 


OP 
dxo 


and finding the smallest positive value of t. For the current approximate solution, 


the characteristics first cross in the region 300 < x9 < 500, 


dx i m(m + 1)AEQ™Brt sin (22g 
_ 5 m+2 


Thus, by setting this equal to zero and re-arranging, we obtain 


Pr m+2 
~200 [D ~ B sin? (=t2y;s00)) 


— 


m(m + 1)AEQ" 2B sin ( 2es-ge0 | 


Hence, we can determine the time the characteristics first cross by finding the 
minimum positive value of t in the region 300 < 2%) < 500. For m = 3, B= Im, 
¢ = 0.4 and Q, = 10, we can determine that the approximate solution is only valid 
until 

t = 238079.124A~' seconds 
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3.6 Numerical Results 


In the next few sections, we use the classic Lax-Wendroff scheme and the flux-limited 
version of Roe’s scheme with the different formulations to obtain a numerical solution 
of either Channel Test Problem A or B. Both numerical schemes use Az = 10m and 
vy = 0.8 and the flux-limited version of Roe’s scheme uses the minmod flux-limiter. 


All schemes use the basic free flow boundary conditions, 


, n+l _ 


= ntl _ won n 
Wag = Woj> Wragg = Wrj Wij 


Wi and Wr; = Wy; 
where i,7 = 1 to 5. Formulation A-CV used a tolerance level of tol = 0.000001 to 
determine when the shallow water equations have reached an equilibrium state. If a 
formulation’s numerical results are not illustrated for a particular scheme then this 


is due to the scheme becoming unstable for this formulation. 


3.6.1 Channel Test Problem A: Numerical Results for a 
Small Bed which is Interacting Slowly with the Water 
Flow 


For the first test case, we simulate a small pulse in the riverbed which is interacting 
slowly with the water flow, where the water flow is moving slowly. To simulate this, 
we use the values A = 0.001, B = 1m and Q = 10. Thus, by using these values 
with the two different schemes, we obtain the numerical results of the classic Lax- 
Wendroff scheme, which are illustrated in Figure 3.7, and the flux-limited version of 


Roe’s scheme, which are illustrated in Figure 3.8. 


From Figure 3.7, we can see that the classic Lax-Wendroff scheme has produced 
poor numerical results for all formulations. The scheme became unstable for 
Formulations A-CV and A-SF and the numerical results of Formulation A-NC were 
completely inaccurate due to spurious oscillations occurring initially in the numerical 


results that were then propagated over the whole domain. This resulted in the 
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B(x,238079s) 


-0.1 
400 500 600 700 


—— Approximate Solution ——A-NC —»— B —&}-C 


Figure 3.7: Comparison of the different formulations for Channel Test Problem A 
using the classic Lax-Wendroff scheme with A = 0.001 and Q = 10 at t = 238079s 


(B). 


B(x,238079s) 


400 500 


A-NC — A-SF x—B —=— C 


Approximate Solution - - - - - - A-CV 


Figure 3.8: Comparison of the different formulations for Channel Test Problem A 
using the flux-limited version of Roe’s scheme with A = 0.001 and @ = 10 at 


t = 238079s (B). 
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total height of the river, the height of the riverbed and the overall velocity 
being greatly reduced (initially, the overall velocity was approximately 1m/s, but 
at t = 238079s, the overall velocity was approximately 0.92m/s). The scheme is 
considerably more accurate with Formulations B and C, but the numerical results 


were still poor due to spurious oscillations being present in the numerical results. 


From Figure 3.8, we can see that the flux-limited version of Roe’s scheme has 
produced accurate numerical results for all formulations except Formulation A-SF 
due to spurious oscillations beginning to occur in the numerical results. For 
Formulation A-CV, the scheme produced numerical results that were very close 
to the approximate solution (for the bathymetry). The scheme produced practically 
identical numerical results for Formulations A-NC, B and C with no spurious 
oscillations present. However, the numerical results of Formulations A-NC, B and 
C were more diffusive than Formulation A-CV due to Formulation A-CV requiring 
considerably less time steps to reach the final time than Formulations A-NC, B and 
C. Formulation A-CV had an overall time step of At ~ 2.9 hours and required 
approximately 3 minutes to iterate the water flow to an equilibrium state each time 
the bed was updated whereas Formulations A-NC, B and C had an overall time step 
of At + 0.7 seconds). Also, notice that Formulations A-NC, B and C have produced 
slightly different numerical results than the approximate solution and Formulation 
A-CV due to the position of the top of the pulse in the bathymetry. Formulations 
A-NC, B and C placed the top of the pulse in the bathymetry slightly to the left of 


Formulation A-CV and the approximate solution. 


To determine if the numerical schemes become unstable for longer computational 
run times, we run the test problem with the same values for longer until ¢ = 150 
hours (540,000 seconds). Thus, we obtain the numerical results of the classic Lax- 
Wendroff scheme, which are illustrated in Figure 3.9, and the flux-limited version of 


Roe’s scheme, which are illustrated in Figure 3.10. 


From Figure 3.9, we can see that the numerical results of the classic Lax-Wendroft 
scheme have become considerably worse for a longer computational run time. All 


variations of Formulation A have now become unstable and Formulations B and C 
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Figure 3.9: Comparison of the different formulations for Channel Test Problem A 
using the classic Lax-Wendroff scheme with A = 0.001 and Q = 10 at t = 150h (B). 
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Figure 3.10: Comparison of the different formulations for Channel Test Problem 
A using the flux-limited version of Roe’s scheme with A = 0.001 and Q = 10 at 


t = 150h (B). 
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have both produced very poor numerical results with spurious oscillations present. 


From Figure 3.10, we can see that the flux-limited version of Roe’s scheme 
has again produced very accurate numerical results for all formulations, with the 
exception of Formulation A-SF that suffered badly from spurious oscillations. Notice 
that the difference of the position of the top of the pulse in the bathymetry between 
Formulation A-CV and Formulations A-NC, B and C has increased. Formulations 
A-CV and A-NC are identical apart from Formulation A-CV iterates the water flow 
to an equilibrium state each time the bed is updated. Thus, the difference is created 
by Formulation A-CV iterating the water flow to an equlibrium state each time 
the bed is updated. If the test problem is run for even longer then the difference 


becomes more significant. 


Hence, the flux-limited version of Roe’s scheme is considerably more accurate 
than the classic Lax-Wendroff scheme for all formulations. Formulation A-SF is 
the least accurate formulation for both schemes and Formulations A-NC and A-CV 
were accurate with the flux-limited version of Roe’s scheme but became unstable 
with the classic Lax-Wendroff scheme. Formulations B and C worked well for both 
schemes but the classic Lax-Wendroff scheme produced spurious oscillations in the 


numerical results. 


3.6.2 Channel Test Problem A: Numerical Results for a 
Small Bed which is Interacting Quickly with the 
Water Flow 


In the previous section we used a value of A = 0.001, which simulated a slow moving 
riverbed. We now use a value of A = 1 to simulate a fast moving riverbed, which 
is now reacting quickly with the water flow. There are a few cases where the bed 
interacts quickly with the water flow, see Perdreau and Cunge [31], but these cases 
are rare. However, it is interesting to study the numerical results of the different 


formulations for this case so that we can determine which formulations can be used 
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with large values of A. Unfortunately when using large values of A, the initial 
conditions used produce an impulsive start for all unsteady approaches due to the 
initial conditions not representing a fast moving riverbed. Formulation A-CV is the 
only formulation that does not produce an impulsive start with the initial conditions 
due to the formulation being a steady approach, which assumes the riverbed is 
interacting slowly with the water flow. At present, obtaining initial conditions for 
a fast moving riverbed can be very difficult. To obtain such initial conditions, 
empirical data could be used but physical cases where the riverbed is interacting 
quickly with the riverbed are very rare. Alternatively, the initial conditions could 
be obtained by iterating an unsteady approach, which includes bed movement, until 
the water flow has settled down, see Hudson & Sweby [21]. 


Now, by using Q = 10, B = 1m and A = 1 with a final time of t = 238s, we 
obtain the numerical results of the classic Lax-Wendroff scheme, which are illustrated 
in Figure 3.11, and the flux-limited version of Roe’s scheme, which are illustrated 
in Figure 3.12. Notice that the pulse in the riverbed has moved at a considerably 


quicker wave speed than with the smaller value of A. 


From Figure 3.11, we can see that by using the classic Lax-Wendroff scheme, 
the scheme became unstable for all of the variations of Formulation A. The scheme 
only remained stable for Formulations B and C, but the numerical results of both 


schemes suffered from spurious oscillations. 


From Figure 3.12, we can see that by using the flux-limited version of Roe’s 
scheme, the scheme became unstable for Formulations A-SF and A-NC. 
Formulations A-CV, B and C have all produced smooth numerical results free 
from spurious oscillations. However, as with the classic Lax-Wendroff scheme, the 
numerical results obtained from Formulation A-CV differed significantly from all the 
other formulations. The pulse in the riverbed has been moved at a slower speed for 
Formulations B and C than Formulation A-CV. The difference between Formulations 
B and C and Formulation A-CV still occurs when the scheme is used with a finer 
mesh, Ax = 0.5, see Figure 3.13. The difference also occurred with the small value 


of A used in the previous section, but the difference was considerably less significant. 
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Figure 3.11: Comparison of the different formulations for Channel Test Problem A 
using the classic Lax-Wendroff scheme with A = 1 and Q = 10 at t = 238s (h & B). 
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Figure 3.12: Comparison of the different formulations for Channel Test Problem A 
using the flux-limited version of Roe’s scheme with A = 1 and Q = 10 at t = 238s 
(h & B). 
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Figure 3.13: Comparison of the different formulations for Channel Test Problem 
A using the flux-limited version of Roe’s scheme on a fine mesh with A = 1, and 
C= 10 abe = 238s (hk & 8): 


The difference is due to the assumptions made when deriving a steady approach. 
The steady approach assumes that the water flow has a negligible effect on the bed, 
which results in the system being decoupled into a water flow approximation, that 
is iterated to an equilibrium state, followed by a bed update. By iterating the water 


flow to an equilibrium state, the equation for conservation of mass (1.1) becomes 


Q,=0 = Q(a, t) =Q. V(z, t), 


where Q, is a constant. Thus, when the water flow has reached an equilibrium state, 
the discharge is constant over the whole domain. Hence, since Formulation A-CV 
iterates the water flow to an equilibrium state, this has the effect of imposing a 


constant discharge. 


Figure 3.14 and Figure 3.15 illustrate the discharge, Q, for each formulation 
compared to the assumed value, Q,. = 10, with A = 0.001 and A = 1 respectively. 


From Figure 3.14, we can see for A = 0.001, the approximate values of the discharge 
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Figure 3.14: Comparison of the discharge for the different formulations with A = 
0.001 at t = 238079s (Q). 
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Figure 3.15: Comparison of the discharge for the different formulations with A = 1 
at t= 2388 (Q). 
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are very close to the assumed values for all formulations. From Figure 3.15 we 
can see for A = 1, the approximate values of the discharge for Formulation A-CV 
are practically identical to the values obtained with A = 0.001 and are close to 
the assumed values. However, Formulations B and C have produced values that 
are beginning to differ from the assumed value. Formulation A-CV assumes the 
water flow is in and equilibrium state, which has the effect of imposing a constant 
discharge. This implies that Formulation A-CV becomes invalid as A — 1 and can 
only be used for small values of A whereas Formulations B and C can be used for 
all values of A. More detailed comparisons of the different formulations with A = 1 
can be found in Hudson & Sweby [21], where it is also illustrated that the numerical 
results of Formulations C and A-CV differ. 


3.6.3. Channel Test Problem A: Numerical Results for a 
Large Bed which is Interacting Slowly with the Water 
Flow 


For the third test problem, we increase the height of the pulse in the bed from 
B = 1m to B = 5m. The riverbed is interacting slowly with the water flow, A = 
0.001, and the discharge is Q = 10. By using these values, we obtain the numerical 
results of the classic Lax-Wendroff scheme, which are illustrated in Figure 3.16, 
and the flux-limited version of Roe’s scheme, which are illustrated in Figure 3.17. 
The numerical results are for t = 5339s and are illustrated with the approximate 


solution, which cannot be calculated beyond this time. 


From Figure 3.16, we can see that the classic Lax-Wendroff scheme has again 
produced spurious oscillations in the numerical results. The scheme has moved the 
pulse at the same wave speed for all formulations, but Formulation A-SF suffers 
the most from spurious oscillations. Formulation A-CV produced results with the 
least spurious oscillations present and the results were similar to the approximate 


solution. 
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Figure 3.16: Comparison of the different formulations for Channel Test Problem A 
using the classic Lax-Wendroff scheme with A = 0.001, B = 5m and Q = 10 at 
t = 5339s (B). 
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Figure 3.17: Comparison of the different formulations for Channel Test Problem A 
using the flux-limited version of Roe’s scheme with A = 0.001, B = 5m and Q = 10 
at t = 5339s (B). 
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Figure 3.18: Comparison of the different formulations for Channel Test Problem 
A using the flux-limited version of Roe’s scheme on a fine mesh with A = 0.001, 
() = 10 end b= bye at t= 53098 (8), 


From Figure 3.17, we can see that the flux-limited version of Roe’s scheme 
produced very accurate results for all formulations, except for Formulation A-SF, 
where the scheme is starting to produce spurious oscillations. The numerical results 
for all formulations were very close to the approximate solution. Figure 3.18 
illustrates the numerical results obtained using the flux-limited version of Roe’s 
scheme on a finer mesh, Az = 0.5, with the different formulations. Here, we can 
see that all of the different formulations with the scheme have produced practically 
identical results. Formulation A-SF with the scheme has produced spurious 
oscillations whereas all of the other formulations produced smooth numerical results. 
Notice that the top of the pulse in the bed has been moved slightly further with 
the different formulations than the approximate solution. This implies that the 


approximate solution has become invalid for B = 5m. 


Hence, the flux-limited version of Roe’s scheme was again considerably more 
accurate than the classic Lax-Wendroff scheme. The numerical results of all 


formulations with the flux-limited version of Roe’s scheme were very similar to the 
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approximate solution if Ax = 10 was used, but the results started to differ if a finer 
mesh was used, i.e. Av = 0.5. Formulation A-SF with both schemes produced the 


least accurate numerical results due to spurious oscillations being present. 


3.6.4 Channel Test Problem A: Numerical Results for a 
Large Velocity with a Small Bed which is Interacting 
Slowly with the Water Flow 


So far, all variations of Channel Test Problem A have been with a small velocity 
resulting in a small Froude number. However, for this test case, we use Q = 50 
with A = 0.001 and B = 1m to simulate a small pulse in a riverbed that is still 
interacting slowly with the water flow, but the speed of the water flow has been 
increased. For this test problem, the Froude number is approximately 0.5, and 
the flow is still subcritical. By using these values, we obtain the numerical results 
of the classic Lax-Wendroff scheme, which are illustrated in Figure 3.19, and the 


flux-limited version of Roe’s scheme, which are illustrated in Figure 3.20. 


From Figure 3.19, we can see that the classic Lax-Wendroff scheme has again 
produced spurious oscillations in the numerical results. For Formulation A-CV, the 
scheme has become unstable and the scheme has produced poor numerical results 
for all formulations, with the exception of Formulation A-NC. Surprisingly, the 
numerical results of Formulation A-NC do not have any spurious oscillations present 


even though the classic Lax-Wendroff scheme is used. 


From Figure 3.20, we can see that the flux-limited version of Roe’s scheme 
produced very accurate results for Formulations B and C but produced spurious 
oscillations for Formulation A-SF. The scheme became unstable for Formulation 
A-CV and with Formulation A-NC, and a small kink has appeared at the top of 
the pulse in the bed. Formulations A-CV and A-NC require an approximation of 
the wave speed of the bed-updating equation, see Section 3.4. The wave speed 


approximation of Chesher et al. (3.14) was used for Formulations A-CV and A-NC, 
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Figure 3.19: Comparison of the different formulations for Channel Test Problem A 
using the classic Lax-Wendroff scheme with A = 0.001 and Q = 50 at t = 1904s 
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Figure 3.20: Comparison of the different formulations for Channel Test Problem 
A using the flux-limited version of Roe’s scheme with A = 0.001 and Q = 50 at 


t = 1904s (B). 
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Figure 3.21: Comparison of the different formulations for Channel Test Problem 
A using the flux-limited version of Roe’s scheme with A = 0.001 and Q = 50 at 
t = 1904s (h+ B). 
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Figure 3.22: Comparison of Formulations A-CV and A-NC for Channel Test Problem 
A using the flux-limited version of Roe’s scheme with A = 0.001 and Q = 50 at 
t = 1904s (B). 
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Figure 3.23: Comparison of the different formulations for Channel Test Problem A 
using the flux-limited version of Roe’s scheme on a fine mesh with A = 0.001 and 
Q = 50 at t = 1904s (h & B). 


but the approximation assumes that the surface elevation is constant, i.e. h © 10—B. 
From Figure 3.21 we can see that the numerical values of the surface elevation 
for the different formulation are not constant. This implies that for Formulations 
A-CV and A-NC, the wave speed approximation of Chesher et al. has become 
inaccurate. In Section 3.4, we also discussed the wave speed approximation of 
De Vries (3.15). Figure 3.22 illustrates the numerical results of the flux-limited 
version of Roe’s scheme with Formulations A-CV and A-NC using either the wave 
speed approximation of Chesher et al. (A-CV/NC-C) or De Vries (A-CV/NC- 
D). The numerical results of Formulations B and C are also illustrated. From 
the results, we can see that Formulation A-NC has produced smoother numerical 
results with the wave speed approximation of De Vries due to the kink at the top 
of the pulse in the bed no longer being present. The scheme became unstable for 
Formulation A-CV with the wave speed approximation of Chesher et al., but with 
the approximation of De Vries, the scheme remained stable and produced smooth 
numerical results. However, the formulation moved the pulse in the bed at a faster 
speed than Formulations A-NC, B and C. This difference still occurs when the 
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Figure 3.24: Comparison of the different formulations for Channel Test Problem 
B using the flux-limited version of Roe’s scheme with A = 0.001 and Q = 10 at 
£= 2508 (2). 


scheme is used with a finer mesh, Ax = 0.5, see Figure 3.23. 


Hence, overall the flux-limited version of Roe’s scheme was again considerably 
more accurate than the classic Lax-Wendroff scheme. Formulations B and C with 
the flux-limited version of Roe’s scheme produced accurate numerical results with 
no spurious oscillations present. Formulation A-NC produced accurate numerical 
results for both schemes and the wave speed approximation of De Vries was the 
most accurate. Formulation A-CV became unstable with the classic Lax-Wendroff 
scheme, and only remained stable with the flux-limited version of Roe’s scheme if 
the wave speed approximation of De Vries was used. Formulation A-SF produced 


completely inaccurate numerical results for both schemes. 
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3.6.5 Channel Test Problem B: Numerical Results for a 
Small Discontinuity in the Bed 


So far, we have only numerically approximated variations of Channel Test Problem 
A. We now use Channel Test Problem B with A = 0.001, Q = 10, Br = 1m 
and Br = 0m to simulate a small sediment bore propagating in a river. By using 
these values, we obtain the numerical results of the the flux-limited version of Roe’s 
scheme, which are illustrated in Figure 3.24, at t = 250h. The numerical results 
obtained with the classic Lax-Wendroff scheme are not illustrated due to the scheme 


becoming unstable for all formulations. 


From Figure 3.24, we can see that the flux-limited version of Roe’s scheme 
produced spurious oscillations in the numerical results for Formulation A-SF. The 
scheme did not produce any spurious oscillations for any of the other formulations. 
Formulations A-NC, B and C produced practically identical numerical results but 
moved the sediment bore slightly further than Formulation A-CV. This is due to 
Formulation A-CV assuming the discharge and surface elevation are constant, which 
is clearly not the case when a discontinuity is present in the riverbed. Notice that 
Formulation B has moved the sediment bore at the same speed as Formulation C 
even for this long computational run time. This implies that Formulation B may 
be written in conservative variable form, but more investigation is required to verify 
this. 


3.7 Summary 


Throughout this chapter, we have discussed five different formulations with two 
different numerical schemes that can be used to numerically approximate the 
equations governing sediment transport. For both test problems, we have seen that 
the classic Lax-Wendroff scheme produced poor numerical results with spurious 


oscillations present. The flux-limited version of Roe’s scheme did not produce 
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spurious oscillations for Formulations B and C, but did with Formulation A-SF. 
Formulation A-CV with the flux-limited version of Roe’s scheme also produced poor 
numerical results for test problems with A = 1 and Q = 50 but produced accurate 
numerical results for small values of A and Q. Formulation A-SF was extremely 
inaccurate due to the numerical results suffering from spurious oscillations if either 
scheme was used. This is due to the formulation having a singular Jacobian matrix. 
If a small value of A is used, Formulations B and C with both schemes suffered 
more from diffusion than Formulation A-CV. This is due to Formulation A-CV 
discretising the water flow and bed update separately resulting in considerably less 
time steps required to reach the final computational run time. Formulation A-CV 
has an overall morphological time step, which is calculated by using 
vAx 3€ Au? 


At = —~ h = 
(e mae where A 7 


Now, 
3€Amax;(|w3}) 
min;(|h]) 


max;,(|A|) = 


and from the numerical results of Channel Test Problem A with Q = 10, B = 1 and 
D = 10, max;(|u|) © 1.1 and min;(|h|) + 9 thus, we obtain 


max;(|A|) = 0.759795A, 


and by using, Ar = 10 and v = 0.8, 


VAL 


ee = 10.52915589A* 
75979BA ~ 1092919989 seconds 


At 


For small values of A, the morphological time step is large, i.e. if A = 0.001 then 
At = 10529 seconds, and as A — 1, At — 10 seconds. This implies that for small 
values of A, Formulation A-CV takes a small amount of time steps to reach the final 
computational run time but as A tends to one, the formulation takes longer due to 
the formulation iterating the water flow to an equilibrium state each time the bed 


is updated. For the other formulations, with A = 0.001 the eigenvalues are 
Ay © —8.28, Ag 0.000772 and A3 % 10.51 
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Figure 3.25: Comparison of the different schemes with Formulation A-NC for 


A 


Channel Test Problem A using A = 1, B= 1 and Q = 10 at t = 238s (B). 


and the eigenvalues associated with the water flow, A; and A3, are the dominant 
eigenvalues and determine the overall time step of the scheme, which is approximately 
At ~ 0.7. As A becomes large, the eigenvalue associated with the bed-updating 
equation becomes the dominant eigenvalue but for realistic test problems, the value 
of A is usually small. Thus a numerical scheme that can be used with a large time 


step is required for all unsteady formulations, such as an implicit scheme. 


Formulation A-NC with the flux-limited version of Roe’s scheme produced 
accurate numerical results for test problems where A is small but as A — 1, spurious 
oscillation started to appear in the numerical results. We can minimise these 
oscillations by using a smaller Courant number, but this can be impractical due 
to long computational run times. As an alternative, we could use the flux-limited 
version of Roe’s Scheme to numerically approximate the shallow water equations 
with 

Brit = BP — s€(ai,.1— 1); (3.17) 
2 2 
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where 
; @+51—vAGh-—a)o? ify 
Gini — 7+ HR Gi — POP if 


T+2 ’ 


> 0 
<0 


to numerically approximate the bed-updating equation for Formulation A-NC. 
Unfortunately, the numerical scheme becomes unstable resulting in the scheme being 


considerably worse than the normal approach (3.12). However, by using a staggered 


approach, 
: (ee vit) Geet —q")o? if vine > 0 (3.18) 
i = n+ n n+ n+ no: n+ 2 ‘ 
“2 a — nl oF vit) Gent —qi"")O? if Vint <0 


we can minimise these spurious oscillations without having to reduce the Courant 
number, see Figure 3.25. Here, we can see that all oscillations have been 
eliminated resulting in a smooth numerical solution, that is considerably more 


accurate than the normal approach. 


Formulations B and C with the flux-limited version of Roe’s scheme started to 
produce different numerical results than Formulation A-CV as A — 1 due to the 
formulation assuming that the discharge and total height of the river are constant 
throughout the domain. Thus, Formulation A-CV should only be used for test 


problems where the riverbed is interacting slowly with the water flow. 


Hence, the flux-limited version of Roe’s scheme was considerably more accurate 
for all formulations. Formulation A-NC produced accurate numerical results for 
test problems where the riverbed interacts slowly with the water flow but became 
unstable when the riverbed interacts quickly with the water flow, unless the staggered 
scheme was used. Formulation A-SF was the worst formulation and both schemes 
with Formulation A-CV became unstable when the riverbed interacts quickly with 
the water flow or the velocity is large. Formulation A-CV produced different 
numerical results to Formulations A-NC, B and C when the riverbed is interacting 
quickly with the water flow and both schemes become unstable when the velocity 
is large. Thus, Formulation A-CV should only be used for test problems where Q 


and A are small whereas Formulations B and C can be used for any test problem 
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with any value of A and Q. In the next chapter, we extend the numerical schemes 


discussed in this chapter and the previous chapter to two dimensions. 
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Chapter 4 


Numerical Schemes for Systems of 
Conservation Laws in 'T'wo 


Dimensions 


So far, we have only discussed sediment transport in one dimension. In this chapter 
we discuss a variety of numerical techniques that can be used to approximate 
the equations governing sediment transport in two dimensions. The classic Lax- 
Wendroff scheme in 2D with either a pointwise or second order accurate source term 
approximation, a flux-limited 2D Lax-Wendroff scheme, a flux-limited 2D version of 
Roe’s scheme and dimensional splitting will be discussed and the numerical results 
compared using a 2D wave propagation test problem. The source terms will be 
constructed to satisfy a 2D version of the C-property of Bermtidez & Vazquez [1]. 
All numerical schemes will be derived so that they can be used to approximate 


systems of conservation laws with source term in 2D, 


Ow  OF(w) _ 0OG(w) 
at = Ow —t—<“‘iéiY 


=H. (4.1) 
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where F(w) and G(w) are the flux-functions and R is the source term, which can 
be re-written as 


R=f+g, (4.2) 


where f and g contains the x or y derivative terms only of the source term 
respectively. To help illustrate the numerical techniques, we use the 2D shallow 


water equations, 


h uh vh 0 
uh | + | hu?+sgh? | + huv = | —ghB, |, (4.3) 
vh : huv : hv? + sgh? ; —ghB, 


where h(x,y,t) is the height of the water above the bottom of the channel (m), 
B(«,y,t) is the height of the riverbed (m) and u(z,y,t) and v(x,y,t) are the 
velocities in the x and y direction, respectively (m/s). Some of the numerical 
schemes discussed in this chapter require the Jacobian matrices of the fluxes, which 


for the 2D shallow water equations are 


oF 0 1 O ye 0 0 
A= =] @=y? 2u 0 ud. B= = —uv , 
Ow Ow 
—w vu c= 0 Qu 


where c = gh, whose corresponding eigenvalues for A are 
Ma=u-c, M=u and MW=ute, 


and for B are 


M=u-—c¢ AS=v and AS H=v+e. 


The eigenvectors are, for A 
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Figure 4.1: Illustration of the initial bathymetry for 2D Test Problem B (B). 


and for B 
0 1 
_ = G — 
= u , ey =| -c and e; = u 
U-—C 0 v+e 


4.1 2D Test Problem: Wave Propagation Test 
Problem 


We use a two dimensional wave propagation test problem, which was discussed 
by LeVeque [26] and Hubbard & Garcia-Navarro [18], for the 2D shallow water 
equations to illustrate the accuracy of the numerical schemes discussed in this 


chapter. For this test problem, the riverbed is fixed. The initial conditions consist 


113 


of 
1.01—B(z,y) if01<2<02,0<y<1 


1— Biz, y) otherwise 


d 


A(x, y,0) = 


ueg,0=0, wie,y,0)—0 


with bathymetry 
B(x, y) = 0.5exp (—50 ((a — 0.5)? + (y — 0.5)”)) , 


which is illustrated in Figure 4.1. We use a gravitational constant of g = 1, which 
was used by LeVeque [26]. For this test problem, the small disturbance in the river 
splits into two waves which propagate in opposite directions. One of the waves 
propagates over the hemisphere that is present in the riverbed. This can cause 
difficulties in obtaining an accurate numerical approximation, depending on how 


the source term is approximated. 


4.2 Conservative Numerical Schemes in 2D 


For 2D systems, there are numerous techniques that can be used to approximate 
the system (4.1). In this thesis, we consider two techniques: a basic finite difference 
scheme and dimensional splitting. Both approaches are discretised on a cartesian 
mesh, see Figure 4.2, where Ax and Ay are the spatial step sizes in the x and y 


* 


direction respectively and F 


hy and Gi, 41 are the numerical flux-function in the 
x and y direction respectively. As with the 1D case, we non-dimensionalise the 
variables, ; A é 
iG 
Pop wap Cap Bap Baz, 
5. ge x. ef 5 a 
Ca Se and wu aa i 
where 


L 
L=maxi(|t1j — £03], |Yi0 — Yiol) and T= 4] . 
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Figure 4.2: The two dimensional mesh. 


denote the non-dimensional coefficients and L is the maximum length of the domain. 
Non-dimensionalising the variables in this way results in the spatial and time step- 
sizes being less than one, i.e. Ax, Ay < 1 and At < 1, which ensures the error 
of a numerical scheme does not grow. For all schemes, we use the basic free flow 


boundary conditions, 


n n+1 __ n 
Wi0 and Wi = Win 


a; 


—_ n+1 __ n n+1 __ 
Wig = Wag. Wriig — Wry ij — 


where 2, 7 = 1 to 5. 


4.2.1 Basic Finite Difference Scheme 


The most basic approach for approximating (4.1) is to construct a finite difference 


scheme of the form 


Wile. ( *, —F ) ~s, (Gris = Cs) + AtR* (4.4) 


iJ tJ i+$5J i-$J ay? 
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where 


At _ At 
~ Ag? 78 Ay 


The numerical flux-functions F* aig and Gi, i and source term approximation R;; 
oy 3 3 


can be obtained by either discretiaing the whole system (4.1) or by discretising 
the system separately in the two coordinate directions. If the system is discretised 
separately in the two coordinate directions, then any of the one dimensional 
numerical fluxes discussed in Chapter 2 can be used and the source term 


approximation is re-written as 


where f;, and g;, are the source term approximations that contain the x and y 
direction components only. To ensure the basic finite difference scheme remains 


stable, we use a variable time step 


vy min(Az, Ay) 
max;,j(|A"|, |A®])’ 


t= 


where A* and A® are the eigenvalues of the Jacobian matrices A and B respectively 
and v is the required Courant number. Unless stated, all versions of the basic finite 


difference scheme discussed in this chapter are stable for v < S. 


4.2.2 Dimensional Splitting Scheme 


Strang [39] derived a dimensional splitting approach that can be used to approximate 
two dimensional systems of conservation laws (4.1). The approach discretises the 


system in parts, effectively solving 


low | OF(w) _ 
20t Ox = oa) 
= Law  dG(w) 
‘WwW WwW 
> Ok =F by =g (4.7) 
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by re-writing the source term as (4.2). The most basic form of dimensional splitting 
is to discretise in two steps, where an approximation of (4.6) is obtained using a half 
a time step, 4!. The results are then used as the initial conditions for (4.7), which is 
then approximated using a half a time step to obtain results at the next time level 


t”*! Hence, we obtain the first order dimensional splitting scheme, 


1 - : : 
wi) = Wij — Sy Ge F. a ) + Saf; (4.8a) 
and 
n 1 at 1 1 
where F* 


ih and Gi, 42 can be any of the numerical flux functions discussed in 
Chapter 2. Unfortunately, the scheme is only first order accurate in time. However, 
Strang [39] extended the approach to be second order accurate in time by discretising 
in three steps, where an approximation of (4.6) is obtained using a quarter a time 
step, 4t. The results are then used as the initial conditions for (4.7), which is then 
approximated using a half a time step, At These results are then used as the initial 
conditions for (4.6) using a quarter a time step to obtain results at the next time 


level t’+!. Hence we obtain a second order accurate dimensional splitting scheme, 


(1) = n Sa * * Sx ok 
Wig ~ Wig ~ > ( +3 Ey + 9 far (4.9a) 
(2) (1) (1) (1) (1) 
Wii = Wij — Sy car = G2) + Sy8; (4.9b) 
and 
nti _ (2) $e (p22) _ pl) Sx o(2) 
Wig = Wig — 5 Get Es) + 5 Lj: (4.9c) 


One advantage of using dimensional splitting is that the CFL condition of the scheme 
is less restrictive than the basic finite difference scheme. However, both versions of 
the dimensional splitting schemes require considerably more computations than the 
basic finite difference scheme. The first and second order versions of the dimensional 
splitting scheme require two and three times as many calculations than the basic 


finite difference approach respectively. To ensure the dimensional splitting scheme 
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remains stable, we use a variable time step 


y min(Az, Ay) 


Ai= 
max;,(|A"|, |A°]) 


Unless stated, all versions of the dimensional splitting scheme discussed in this 


chapter are stable for v < 1. 


4.2.3. 2D C-Property 


The C-property discussed in Section 2.2.2 can be extended to two dimensions. 


Consider the 2D shallow water equations for the quiescent flow case, 
u(z,y,t)=0, v(a,y,t)=0 and A(a,y,t)=D-B V(z,y,t). 
For this stationary case, w; = 0 thus, the flux functions and source terms balance 
F(w), + G(w), =R. 


In the 1D case, we only needed to balance two terms but here, there are now three 
terms that need to balance, which can be considerably harder to obtain numerically. 
However, by re-writing the source term as (4.2) and splitting the equation in two, 
we obtain 

F(w), =f and G(w),=g. 


Here, the fluxes must balance with the source term approximations. To obtain this 
condition numerically, a scheme must balance the numerical fluxes with the source 


term approximations, 


i+3,j _ a ij+s — wjJ-s 


If the source term approximations balance with the numerical fluxes, then the 


numerical scheme satisfies: 
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e the approximate C-property, if the numerical scheme is accurate to the order 


O(Az?), when applied to the quiescent flow case; 


e the exact C-property, if the numerical scheme is exact when applied to the 


quiescent flow case. 


If a numerical scheme does not satisfy the C-property (exact or approximate) then 


spurious waves may occur in the numerical results, as illustrated in one dimension. 


4.3. Lax-Wendroff Scheme in 2D 


The classic 2D Lax-Wendroff scheme is a basic finite difference scheme (4.4), where 


the whole system is discretised, and is derived by using the Taylor’s series expansion, 


2 


At 
wott a y] w; ai At(we);; + (wu)? + O(At?) 


4,7 ) 


to obtain the numerical flux-functions 


* 1 nm mr n 
Piya y ~ 9 (Eeag ha) SA (3s 


i+3,j 


= + (Aha (aya — GR j-1) + AP; (Gar — G21) (4.10a) 


_F",) 


i+1,7 


(Giiai + GR) — 5 BM (Gia — GRy) 


_= (Beja (FRager — Feige) + BR (Feng —Fii,;)). (4-10) 


The scheme requires an approximation of the Jacobian matrices, which can be 


approximated by averaging the neighbouring cells, 


ie. net ee hs eS Gh ue 
AM a (Huy and B?,,,=B( 4), 
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For the source term approximation, we use a pointwise approach 


0 
Rij = | ake (hag + hRay) (BRag — BEay) (4.11) 
—ghy (AB ja + AB )-1) (BR — BE-1) 


but as with the 1D case, this is a crude approximation of the source term since it 
does not make the scheme satisfy the C-property. The classic 2D Lax-Wendroff is 


i 
stable for v < Fe 


Figure 4.3 and Figure 4.4 show the numerical results obtained using the classic 
2D Lax-Wendroff scheme (4.10) with the pointwise source term approximation (4.11) 
at t = 0.1 for the 2D Test Problem. The scheme was used with Ax = Ay = 0.001 
and v = 0.4. Here, we can see that even though the test problem has not been run 
for very long, the scheme has produced very poor numerical results with spurious 
oscillations present. Even though the water surface is supposed to be smooth away 
from the disturbances, a hemisphere has appeared on the water surface just above 
the hemisphere in the riverbed. The poor results obtained are due to the pointwise 


source term approximation which results in the scheme not satisfying the C-property. 


In 1D, we derived a second order accurate source term approximation for the 
classic 2D Lax-Wendroff scheme, which made the scheme satisfy the C-property. A 
2D version of this scheme can be obtained by including the source term in the Taylor 
series expansion when deriving the classic 2D Lax-Wendroff scheme to obtain 

Ae 


R’ = AiR + = (R, ~ (AR), = (BR), ) 


Now, by re-writing the source term using (4.2), we obtain 


ri 


R’ = At( +2) - 5 ((an), + (Be),) - 


—— (Ag). + (BE),) 


where the term AFR, has been omitted for convenience and is discretised separately. 
Even though this term has been omitted, the source term approximation still ensures 


the scheme satisfy the C-property. Hence, we obtain a second order accurate source 
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Figure 4.3: Numerical results of the classic 2D Lax-Wendroff scheme at t = 0.1 
(A+ B). 
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Figure 4.4: Numerical results of the classic 2D Lax-Wendroff scheme at t = 0.1 
(h+ B). 
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term approximation, 


‘ 1 n me iy 
R;; = Ag ((1 = SpAj1 a) fi J + (I + s2At 1 ;) ae 


’ nm 
a 2Ay ((1- a) Bi j42 a (1 + = :) g", 3) 


— Sry (As )iiay — (Agliay) — axe (BN - BAA): 


which can be re-written in a more convenient way as 


i _ oF = =F 
AR, = 5. (P52, +f1,) +4 (era, teh)» (4.12) 
where i 
a n Sy n n 
fa j 9 (1 a SA je ;) fia, 9 ((Agyia, + (Ag); i) 
and 


ching (eBay) aha 7 (Bd HN). 


For the 2D shallow water equations, we obtain 


0 
fy = = (hea - Ks ) (Bhig - Br) ? 
0 
0 
= — 5h; (BRag Bris) ? 
0 
Si 54) 
<5 (he jas + hi ‘;) (BP ij+1 Br) 
and 
Bij = 
—$h2; (BRiu1 — BBs) 
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The AER, term, which was omitted earlier, can be included by using the chain 


rule R; = Rww;. Thus, we obtain a semi-implicit version of the 2D Lax-Wendroft 


scheme 
At ~ 
n+l n n * * 
a i Wij ee ci 7 a yy] (Fas Fi) 
eee lee ee , es 
= hi - = (Rw); i (G; Naa Gh 1) + At hi - = (Rw) | Re. 


where for the 2D shallow water equations 


0 0 0 
(Rw) i; a "TAs (Brag - Be i) 0 0 
= 5ag (Bria — Bp) 0 0 


To illustrate the accuracy of the second order accurate source term approximation, 
we use the 2D Test Problem. Figure 4.5 and Figure 4.6 show the numerical results 
obtained using the classic 2D Lax-Wendroff scheme (4.10) with the second order 
accurate source term approximation (4.12) at t = 0.7 for the 2D Test Problem. 
The scheme was used with Ax = Ay = 0.001 and v = 0.4. Here, we can see that 
even though the numerical results are considerably more accurate than using the 
pointwise source term approximation (4.11), spurious oscillations have started to 


appear in the numerical results. 


4.4 Flux-Limited Lax-Wendroff Scheme in 2D 


The Lax-Wendroff scheme suffers from spurious oscillations, which can be minimised 
by adapting the scheme so that it satisfies the TVD property. Flux-limiter methods 
are used to obtain a high resolution scheme by constructing numerical flux functions 
of the form 


TVD FO SO FO TVD FO FO 
FIVE, = FFG + @ (FS, - FPG.) and GTY2 = GPO +6 (GS, - GFo ,| 
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Figure 4.5: Numerical results of the classic 2D Lax-Wendroff scheme with the second 
order accurate source term approximation at t = 0.7 (h + B). 
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Figure 4.6: Numerical results of the classic 2D Lax-Wendroff scheme with the second 
order accurate source term approximation at t = 0.7 (h + B). 
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where F°?,. and G°? , are second order numerical fluxes, F*% . and G¥° , are 
i455 Lj+s +55 Lj+s 

first order numerical fluxes, ® = diag(®;,) and ©, is a flux-limiter. By using the 

classic Lax-Wendroff numerical flux (4.10) with the first order upwind numerical 


flux, 
ay ji 9 (Ey + FE; i) _ 5lAia (why _ wij) y 


we obtain the flux-limited Lax-Wendroff numerical flux 


1 i n 
Ft, (F,,,+F",) - 5 (X|A"|LX™) 7, (wi; — wz,) 


i+h gj — 9 
S 
— Pia y (Atay (Gia s41 — Giag-1) — Aly (Ghizt — Giza), (4-188) 


where X is a matrix containing the right eigenvectors, ef, of the Jacobian matrix 


associated with F, AY = diag(Af) is the diagonal matrix of eigenvalues, Af, and 
L = diag (1 — ®(6,)(1 — |X|) . 


Similarly, for G we obtain 


G; = 5 (Ghia + G? “ = 


atk (Y|AC|LY- ) aad (wear — Wis) 


1 
2 
a 7 vi i i 7 
By 54 1] (Fe t1,j7+1 na F; gay) — Bi; (Faiz — i. i) (4.13b) 


a Be aitd ( 


where Y is a matrix containing the right eigenvectors, e@, of the Jacobian matrix 


associated with G, AG = diag(A@) is the diagonal matrix of eigenvalues, \@, and 
L = diag (1 — ®(0,)(1 — |A¥7])) . 


We require a source term approximation for the numerical scheme that makes the 
scheme satisfy the C-property and adopt an approach discussed by Hubbard & 
Garcia-Navarro [18, 19], where the flux functions and source term approximations 
are balanced so that the numerical scheme satisfies the C-property. As with the 1D 
case, we construct source term approximations of the form (4.5) where 


_ TVD 4 TVD TVD _ ¢FO n 
fi; al EW Patel ? fag fag Pag 


Ga ~ fF oF 


i+53 i+4 
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ee ave 


= TVD 
ij = 8; ijt+s + 8i5-4 


and Brit = Bij + Py Cae — Bij ayy 

where the superscripts FO and SO denote first and second order accurate source 
term approximations respectively. In the previous section, we derived a second 
order accurate source term approximation (4.12) that makes the scheme satisfy 
the C-property. We use this approximation with the first order approximation of 


Bermtidez & Vazquez [1] where 


£4 = 5 (@FIAIA) #): 


fh 
I 5d 


to obtain a flux-limited second order approximation of the source term 


Pf, of; 


oJ i+3,j i-4,j? 


+ = if F\-1 F -1 e Sy n n 
fia g= 5 (KX (T+ (4") AFL) X Bas = BP (Ag) + (AB) - 
Similarly, for g we obtain 


gi; = Bi j4a t gy ig}? 


where 


(v (1 = (Ae) |AC|L) Y's) Fa a re (BEY. + (Bf); ) 


ro 
G 
+ 
I 
Nol rR 


To illustrate the accuracy of the flux-limited 2D Lax-Wendroff scheme, we use the 
2D Test Problem. Figure 4.7 and Figure 4.8 show the numerical results 
obtained using the flux-limited 2D Lax-Wendroff scheme at t = 0.7 for the 2D 
Test Problem. The scheme is used with spatial step-sizes Ar, Ay = 0.01, v = 0.4 
and the minmod flux-limiter. Here, we can see that the scheme has produced very 
smooth numerical results, which are in agreement with the results found in Hubbard 
& Garcia-Navarro [18]. The scheme is considerably more accurate than the classic 


Lax-Wendroff scheme with either source term approximation. 
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Figure 4.7: Numerical results of the flux-limited 2D Lax-Wendroff scheme at t = 0.7 
(A+ B). 
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Figure 4.8: Numerical results of the flux-limited 2D Lax-Wendroff scheme at t = 0.7 
(A+ B). 
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4.5 Roe’s Scheme in 2D 


In two dimensions, the flux-limited 2D version of Roe’s scheme is obtained by 
discretising the system separately in the two coordinate directions. Hence, we obtain 


the flux-limited version of Roe’s 2D numerical fluxes 


P 


* u n mr 1 ay \ 
Pia = 5a +P) —5 > [af ARG - Dawe ek]... (414a) 
- 9 5 = +5, 
and 
* l n n 1 : AG|\G G 
Gijng = 3(Glin + GE) — 5 > [AFAR — 9062) — bE DDE,» a) 
where PF) 
On JI4h 3 
Vp = SXk, 6;. = (af) _ [= i — sgn(v}, ig jt 
kith,j 
(ax); J+4 
0G = 7 - =, J = 7 —sgnlv )i ) 
‘ (OR )ij43 ieee 


and ®; can be any of the flux-limiters listed in Table 2.1. The values of Mg and &; 
are the eigenvalues and eigenvectors of the Roe averaged Jacobian matrix, A or B, 
and a, are the wave strengths, which are determined from the decomposition, see 
Section 2.7 for more details. The superscripts F’ and G denote if the terms are for 
Jacobian matrix A or B, respectively. Glaister [11] determined the following Roe 


averaged Jacobian matrices for the 2D shallow water equations, 


0 1 0 0 0 1 
A=| @-# 2% 0 and B=] —-a o @ 
it 80 tt 2 0 2 


where ¢ = 4/ gh, whose corresponding eigenvalues for A are 


W=G-6 M=% and MW =%4+6 
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and for B are 
MS =G-& AMG =H and M=H+E 


The eigenvectors are, for A 


1 0 1 
a =| feel, = =| 0 and é6f =] @+é |, 
Uv c Uv 
and for B 
1 0 1 
= tl , C=) = and ef = u 
oC 0 We 


The wave strengths are for A 
il pl e 
Ge sah + = (uAh — A(uh)) and @, = z (A(vh) — oAh), 
where Aw = wr, — Wz,;, and for B 


Gt3=5 sAh + 5 = (BAA — A(vh)) and 4% = 


where Aw = Wi p — Wiz. The Roe averages are 


VhretVht | Vhrt+ Vat 


where R and L denote the right and left components in either the direction x or y. 


foot] 


C= 


~ VA vh Vvh vh ~ 1 

RUR + LUL RUR+ LUL sae OS : (ia + he), 
For the source term approximation, we also discretise separately (4.2) to obtain 
f= Tag tty One By —2 at Es (4.15) 


its J 1,j— 5 


where 


C. “3 [AP EF (1 + san(AZ)(1 — @e(1 — EI) 


i+h,j 
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n i ee 7 
ef = 5 > [Bek + sen AG) - & C1 - WE) 
k=1 


ijt 


and the values of 6, are the coefficients of the decomposition of the source terms 

onto the eigenvectors of the characteristic decomposition, see Section 2.7 for more 

details. For the 2D shallow water equations we obtain 
~ CAB 


A= a” Bz =0 and —— 


cAB 
2 


for both F and G. By using this source term approximation, the scheme satisfies the 
C-property. The flux-limited version of Roe’s 2D numerical fluxes (4.14) with source 
term approximations (4.15) can be used with the basic finite difference scheme (4.4) 
and the first and second order accurate in time dimensional splitting schemes (4.8) 
and (4.9) respectively. In this thesis, we denote the flux-limited version of Roe’s 2D 
numerical fluxes with the basic scheme as the B-FLR scheme and the second order 


dimensional splitting schemes as the DS2-FLR scheme. 


To illustrate the accuracy of the basic scheme (B-FLR) and dimensional splitting 
scheme (DS2-FLR) we use the 2D Test Problem. The schemes are used with 
Az, Ay = 0.01, v = 0.4 and the minmod flux-limiter. Figure 4.9 and Figure 4.10 
show the numerical results obtained using the basic scheme (B-FLR) at t = 0.7. 
Here, we can see that the B-FLR scheme has produced numerical results similar to 
the flux-limited 2D Lax-Wendroff scheme, see Figure 4.7 and Figure 4.8, but the 
numerical results are slightly more accurate. The results are also similar to those 
obtained in Hubbard & Garcia-Navarro [18]. 


Figure 4.11 and Figure 4.12 show the numerical results obtained using the 
dimensional splitting (DS2-FLR) scheme at t = 0.7. Here, we can see that the 
scheme has produced slightly better numerical results than using the basic finite 
difference scheme. The results are also similar to those obtained in Hubbard & 


Garcia-Navarro [18]. 
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Figure 4.10: Numerical results of the basic scheme (B-FLR) at t = 0.7 (h + B). 
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Figure 4.11: Numerical results of the dimensional splitting scheme (DS2-FLR) at 
t=0.7 (h+ B). 
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Figure 4.12: Numerical results of the dimensional splitting scheme (DS2-FLR) at 
t=0.7 (A+ B). 
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4.6 Summary 


Throughout this chapter, we have discussed a variety of numerical techniques that 
can be used to approximate the equations governing sediment transport. We used 
the 2D shallow water equations to help derive and illustrate the accuracy of the 
numerical schemes. We showed that the classic 2D Lax-Wendroff scheme suffered 
badly from spurious oscillations, which were considerably worse when a pointwise 
source term approximation was used. The second order accurate source term 
approximation eliminated some of the spurious oscillations and greatly improved the 
accuracy of the scheme, but spurious oscillations were still present. A flux-limited 
2D Lax-Wendroff scheme was derived that was very accurate with no spurious 
oscillations present in the numerical results. The flux-limited version of Roe’s 
scheme was adapted to two dimensions and the basic finite difference scheme and 
dimensional splitting schemes were discussed. The numerical results obtained using 
the flux-limited 2D version of Roe’s numerical fluxes were the most accurate with 
the second order accurate dimensional scheme producing the best set of numerical 


results. 


In the next chapter, we adapt the different formulations discussed in Chapter 3 
to two dimensions and use the basic scheme (B-FLR), the second order accurate 
dimensional splitting scheme (DS2-FLR) and the 2D Lax-Wendroff scheme with 


pointwise source term approximation to obtain a numerical approximation. 


133 


Chapter 5 


Numerical Formulations for 
Approximating the Equations 
Governing Sediment Transport in 


Two Dimensions 


In the previous chapter we discussed a variety of schemes that can be used to 
approximate systems of conservation laws in two dimensions. In this chapter, we 
discuss how to approximate accurately the equations governing sediment transport 


in two dimensions, which comprise of the equation for conservation of mass, 


Oh | O(uh) 
Ot Ox 


_ O(vh) 
ay =, (5.0) 


the equation for conservation of momentum in the x-direction, 


O(uh) 
Ot 


O (hu? + igh? O(huv 
( oe ) | mf se (5.2) 
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the equation for conservation of momentum in the y-direction, 


O(vh) zr 0 (huv) ‘. 0 (hv? + $gh?) 


ey Da Oy ==9iB, (5:3) 
and the bed-updating equation, 
OB An) | Aga) _ 
a g an & op = 0, (5.4) 


where € = = and ¢ is the porosity of the riverbed. Here h(x,y,t) denotes the 
height of the water above the bottom of the channel, B(x,y,t) denotes the height 
of the riverbed, u(x, y,t) and v(x,y,t) are the velocities in the x and y direction 
respectively and qi(u,v,h) and qo(u,v,h) are the total (suspended and bedload) 
volumetric sediment transport rates in the x and y direction respectively. The 
sediment transport fluxes are not direct functions of B and are more complex in two 


dimensions. 


In Chapter 3, the steady and unsteady approaches were discussed and used 
to derived five different formulations. Here, we extend the two most accurate 
formulations to two dimensions. ‘To obtain an approximation of the different 
formulations some of the different schemes discussed in Chapter 4 are adapted. 
In this chapter, the classic 2D Lax-Wendroff scheme (4.10) and the basic (B-FLR) 
and dimensional splitting (DS2-FLR) schemes, which both use the flux-limited 2D 
version of Roe’s numerical flux (4.14), are adapted. A 2D Channel Test Problem 
is then used to illustrate the accuracy of the two formulations and the different 


numerical schemes. 


5.1 Different Formulations 


In Chapter 3, we derived five different formulations that can be used to approximate 
the equations governing sediment transport in one dimension. In this section, we 


extend some of the formulations to two dimensions: 
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e Formulation A-CV will be extended to 2D even though the numerical results 
sometimes differed from the other formulations in 1D. This is the formulation 
currently used in industry and is based on the steady approach, where the bed 


is assumed to have a negligible effect on the water flow. 


e Formulation A-NC produced poor numerical results for some of the test 


problems in 1D and therefore will not be extended to 2D. 


e Formulation A-SF produced poor numerical results with spurious oscillations 


present for all schemes in 1D and will not be extended to 2D. 


e Formulation B cannot be extended to 2D since the equation for conservation 
of momentum in the x direction (5.2) can be re-written in non-conservative 


variable form as 


u Ea + g(h4 a) + vuy = 0, 


x 
but cannot be written in conservative form due to the vu, term. This is also 


the case with the equation for conservation of momentum in the y direction. 


e Formulation C will be extended to 2D and was one of the most accurate 
formulations in 1D. This formulation is based on the unsteady approach, where 


the water flow and bed update are discretised simultaneously. 


Thus, Formulation A-CV, which is a steady approach, and Formulation C, which is 
an unsteady approach, will be extended to two dimensions. The two dimensional 


sediment transport fluxes must first be derived for the formulations. 


5.1.1 Sediment Transport in 2D 


In two dimensions, the sediment transport flux is treated as a vector, it’s components 


being given by 
_ vq 
( aa) )= [U]’ 
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where the magnitude of the velocity is 
|U| = Vu? + v? 


and |q| is the magnitude of the sediment transport flux, which is determined from 
any of the one dimensional sediment transport fluxes discussed in Section 1.4. For 


example, if the sediment transport flux (3.3) is used in two dimensions then 


aT 
2 


lq) = AJU|" =A (u? + v’) (5.5) 
Hence, we obtain the sediment transport fluxes in the x and y direction 
qi(u,v) = Au (u? + y2) 20-0) and qo(u,v) = Av (u? =" y2) 2m) (5.6) 


respectively. Notice that in 2D, the sediment transport fluxes have become more 
complicated, which leads to even more difficulties in obtaining an accurate numerical 
approximation. In this thesis, all formulations discussed in 2D will be based on the 


sediment transport fluxes (5.6) with m = 3, 


gi(u,v) = Au(u? +?) and q(u,v) = Av(u? + v”). (5.7) 


5.1.2 Formulation A 


Even though the numerical results obtained from Formulation A-CV differed from 
the other formulations for certain test problems in 1D, we extend the formulation 
to 2D as this is the approach currently used in industry. Formulation A is based 
on the steady approach, which assumes that the bed has a negligible effect on the 
water flow and decouples the system. This formulation was pioneered by Cunge et 


al. [3] and consists of 
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e a “fixed-bottom step”, where the 2D shallow water equations, 


h uh vh 0 
uh | +} hu?+ sgh? -- huv =| —ghB, |, 
vh |, huv : hv? + sgh? ; —ghB, 


are iterated to an equilibrium state whilst keeping the riverbed fixed. 


e a “changing bottom step”, where the riverbed is updated whilst keeping all 


other variables fixes. 


The 2D shallow water equations are iterated to an equilibrium state each time the 
bed is updated and the overall time step is the morphological time step of the bed- 
updating equation. Numerically, an equilibrium state in two dimensions is obtained 
when 
lw —wy,|<tol V(i, 7). 

In two dimensions, this condition is hard to obtain, especially for low tolerance levels 
as there are a considerable amount of points to satisfy. Thus, the numerical schemes 
have a cut off time, where if the tolerance level hasn’t been obtained but the water 
surface is calm, the water flow is assumed to have reached an equilibrium state. This 


cut off time varies considerably depending on which test problem is used. 


For Formulation A, the Jacobian matrices of the 2D shallow water equations will 
be required and are discussed in Chapter 4. An approximation of the wave speeds 


of the bed-updating equation are required, 


v6] ana a = 6 [20], 


OB 


which can be difficult to obtain. In this thesis, we extend the analytical approach, 
which was discussed in Section 3.4, to two dimensions. The analytical approach 
assumes that the discharge and the surface elevation are constant, i.e. 

Qu 


A(z, yt) =D = BG, 1)) Ola.) = h(x, y,t) 


Qu 


and v(2,y, 2) = h(x, y,t)’ 
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and then re-writes the velocities in terms of B, 
u(z,y,t)=Qu(D-—B)" and v(z,y,t) =Q,(D—B)", 


so that the sediment transport fluxes can be re-written in terms of B. For example, 
consider the sediment transport fluxes (5.6) which can now be re-written in terms 
of B, 


qu(u,v) = Au(u2 + v2)2"-) = AQ, (Q2 + Q?2) 2") (D— BY-™ 


and 
(u,v) = Av(u? +09)" = AQ, (Q2 + Q2)2"" (D- BY, 


and by differentiating with respect to B, we obtain analytical approximations of the 


wave speeds, 


(u? + v”) am) and yo (u? + v”) gin dy (5.8) 


— A€mv 
OR 


5.1.3. Formulation C 


A two dimensional version of Formulation C, which is an unsteady approach, can 
be obtained by re-writing the equations for conservation of momentum in the x and 


y directions as 
1 
(uh): + iu’ ss sgh Si onb + (huv), = 9Bhy (5.9) 


and 


iL 
(uh), + (huv), + i’ + 5gh + onb = opi, (5.10) 
y 
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respectively by using the chain rule. Now, combining (5.1), (5.9), (5.10) and (5.4) 


into system form, we obtain 


h uh vh 0 
uh ic hu? + gh? + ghB n huv _ gBhy 
vh huv hv? + sgh? + ghB gBh, 
B |, fn . qe , 0 


(5.11) 


If the sediment transport fluxes (5.7) are used with Formulation C, then the 


Jacobian matrix of F is 


0 1 0 0 
OF gih+B)—u? 2u 0 gh 

A =| ——.: ; 
Ow —uv vu O 


—ud — ve de OQ 


where d = 4° (3u2 + v2) and e = 248. One of the eigenvalues of the Jacobian is 
h h 


AY = u and the other three are obtained by solving the cubic 
P(A, w) = »° — 2ud? + [u? — g(h + B+hd)] A+ ghud = 0. 


In Section 3.1.3, we proved that the roots of this polynomial are always real and 
unequal and deduced that we can use (3.7) to determine the values of the roots. 
Once the roots of P(A,w) have been obtained, they are used to determine the 
eigenvectors, 

1 

rx 

v ; 
u? — g(h + B) + (Ag — 2u)Ax 

gh 


—_— 
e, = 
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for k = 1,2,3 wherea#k+#b and ife £0 then 


0 
. 0 
ef = “= ~ (n+ B) otherwise ef = i 
e€ 
1 0 
—-(h+B 
(+B) 
Similarly, for G the Jacobian matrix is 
0 0 1 O 
B= OG _ —uv v u O 
Ow g(h+ B)—v* 0 2 gh 
—vd — ue e d 0O 


where d = a (3u2 + u?) and e = BABU One of the eigenvalues of the Jacobian is 


AG = v and the other three are obtained by solving the cubic 
P(A, w) = A° = 207 + [v? —gh+BH+ hd)| A + ghud = 0. 


Once the roots of P(A,w) have been obtained, they are used to determine the 


eigenvectors, 
1 
u 
ef = Ar , 
v? — g(h + B) + (Ag — 20) Az 
gh 


fork = 1,2,3 wherea#k +b and ife £0 then 


; 0 
(A+B) 
Ws 1 
e = he otherwise ef = 
VU 0 
1 

= (h+B 0 

r(h+B) 
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Now that we have discussed Formulations A and C, in the next section we adapt 
the classic 2D Lax-Wendroff scheme and the 2D variations of the flux-limited version 


of Roe’s scheme. 


5.2 Adaptation of the Classic 2D Lax-Wendroff 


Scheme 


The classic 2D Lax-Wendroff scheme has numerical flux functions (4.10) and the 
source term is approximated using a pointwise approach, 

R; = R(w;';), 
and is commonly used in industry to approximate the equations governing sediment 
transport in two dimensions. As with the 1D case, the classic 2D Lax-Wendroff 
scheme suffers badly from spurious oscillations occurring in the numerical results. 
The pointwise source term approximation results in the scheme not satisfying the 


C-property. 


5.2.1 Formulation A 


Formulation A is a steady approach, where the system is decoupled into a water 
flow approximation followed by a bed update. The classic 2D Lax-Wendroff scheme 
can be adapted for this formulation. In the previous chapter, the 2D shallow water 
equations were approximated using the classic 2D Lax-Wendroff scheme, see 
Section 4.3 for details. However, the scheme needs to be adapted to approximate 


the bed-updating equation (5.4), which is easily obtained 


rete — Bi, — €Sz (iva, _ (a1)i_3;) — Sy (a2) _ (2);;-3) , (5.12) 


142 


where the numerical fluxes are 


* 1 n n Sx n n 
(ina a) (Caer + (11)f;) — ieee (Chee ~ (11) ;) 


s 
_ = Coen ((q2)Praj41 = (q2)Pr1 5-1) aI Ni (Chere = (q2)?j-1)) 
and 
(92)5 44 = 3 (Cee + (q2)";) _ ee ((g2) P44 a (q2)?,) 


_ = ei (Caer - Chee as x (Chee a (q)P-4;)) 


and where AF and A® are the analytical approximations of the wave speeds of the 


bed-updating equation (5.8). 


5.2.2 Formulation C 


The classic 2D Lax-Wendroff scheme (4.10) can be used as written to approximate 
Formulation C. The Jacobian matrices of the formulation are given in Section 5.1.3 


and the source term is approximated with a pointwise approach, 


0 
KT. = 4Ax (BR 23 BP y 5) (Mee hey i) 
re (Bris + Bh -1) (hija nest) 

0 


5.3 Adaptation of the Basic and Dimensional 
Splitting Schemes 


Both the basic scheme (B-FLR) and the second order accurate dimensional splitting 
scheme (DS2-FLR), which both used the flux-limited 2D version of Roe’s numerical 
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fluxes (4.14), produced very accurate numerical results in the previous chapter. This 
is mainly due to both schemes using the source term approximation (4.15), which 
results in both schemes satisfying the C-property. Thus, in this section, we adapt 


both numerical schemes to approximate Formulations A and C. 


5.3.1 Formulation A 


Formulation A is a steady approach, where the system is decoupled into a water 
flow approximation followed by a bed update. The flux-limited 2D version of 
Roe’s scheme can be adapted for this formulation. In the previous chapter, the 2D 
shallow water equations were approximated using both the basic scheme (B-FLR) 
and dimensional splitting scheme (DS2-FLR), see Section 4.5 for details. However, 
the schemes needs to be adapted to approximate the bed-updating equation (5.4). 


Thus, we obtain the following numerical fluxes for the bed-updating equation, 


(anyaay = 2 (a)bay + Qh) 
- ; Ma] 1-9 (825) (- [Aa g])) (Bag — Bh) -13a) 
and 
(92)5 44 = : (Cee + (q2)7;) 
= ; M43] (1- © (66,2) (1 — [eS .a])) (Bt — Bh), 6-18b) 
where Be Be 
Vid g - SeNyt gp Wri4d > Sys O34 - Bi ~ Br, 
ae a Pepe f=it-sgn Faas) J =j —sgn (%,1) 


and ® can be any of the flux-limiters listed in Table 2.1. Here, A” and \© are the 
analytical approximations of the wave speeds of the bed-updating equation (5.8). 


The numerical fluxes of the bed-updated equation (5.13) can be used with both 
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the basic finite difference scheme (4.4) and the dimensional splitting schemes (4.8) 
and (4.9). In this thesis, if the basic scheme (B-FLR) or the dimensional splitting 
scheme (DS2-FLR) are used then the numerical fluxes (5.13) are used with (4.4) or 
(4.9) respectively. 


5.3.2. Formulation C 


Formulation C can be used as written with both the basic scheme (B-FLR) and the 
dimensional splitting scheme (DS2-FLR). Both versions require the Roe averages 
of the system. For the sediment transport fluxes (5.7), the Roe averaged Jacobian 


matrix for F is 


0 1 0 0 

eo gh+B)—i 2% 0 gh 
=UuUU U U 
-~tid-té d é 0 


One of the eigenvalues of the Jacobian is \F = % and the other three are obtained 


by solving the cubic 
PO) = 38 — 202 + i? —g(h + B+hd)| \ + ghad =0. 


The roots of P()) are determined by using the approach discussed in Section 3.1.3. 
Once the Roe averaged eigenvalues have been obtained, they are used to determine 


the eigenvectors, 
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fork = 1,2,3 wherea#k +b and ife £0 then 


1 0 
u 0 
ef = f otherwise é{ = a 
=(h+B 0 
a ) 
e . if. 2 
where f = 0 — = (i + B). The wave strengths are 
he 
e if€ £0 then 
er Vk 
= a Se 
(f — B)(Ag — Aa) (Ak — Av) 
where 
tbe = |(2i — Aa — Aa)iid — (+ B)gi + (Ande + g(h + B) — a) f] An 
+ (f —6) [ghAB + (2% —1y— A)A(uh)] - G—Aa)G-} 


e if€=Othen 


; (a\o + g(h + B) — ti?) Ah + (2% — Aq — Xx) A(uh) + ghAB 
k — mF ~ ~ ~ 

(Az — Aa) (Az — Av) 
for k = 1,2,3, wherea 4k #b and 


For the source term, 


2% — Xa — Ax) GBA , 
= ( . Aq As)9 - and (4, =0, 
(Az — Aa) (Az — Ab) 
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where k=1,2,3 anda4kF b. 


Similarly for G, the Roe averaged Jacobian matrix is 


0 0 1 0 

BE ; ww v wu 0 . 
gh+B)—-t 0 20 gh 
~td — tié éd 0 


One of the eigenvalues of the Roe averaged Jacobian is AG = 6 and the other three 


are obtained by solving the cubic 
P(X) = 33 — 2052 + Gi —g(h+ B+ hd] 1+ ghod =0. 


The roots of P(X) are determined by using the approach discussed in Section 3.1.3. 
Once the Roe averaged eigenvalues have been obtained, they are used to determine 


the eigenvectors, 


1 
Uu 
ey _ Ne } 
6? — g(h+ B) + (Ag — 20) 
gh 


. 0 
Ef = / otherwise éf = ; : 
oe 
a0 + B) 0 
= 


Oe = (fF — @) |ghAB + (26 — Xa — Ay)A(wh)| — (6 — Aa) — No) A(uh) 
oa — (h+ B)git + (Addy + g(h + B) — 6) f| Ah 


+ |(26-X_— 4a) 
for k = 1,2,3, wherea 4k 4b and 
A(uh) — uAh 
ag = >... 
f-wu 


e ifé=0 then 
(arp + g(h + B) — 6?) Ah t+ (26 — Aq — Ax) A(vh) + ghAB 


(Ax = Aa)(Ak = Av) 


ay = 


for k = 1,2,3, wherea #44 4b and 


For the source term, 

> 6—\,—A»)gBA 3 

6 ee ot V0 
(Az — Aa) (Az — As) 


where kK=1,2,3 anda4kF b. 


The Roe averaged values are 

. WVhrurt+Vhiv, . Vhrur+Vhpu, +; 1 

Uu= >» U= 5 h=-=(hre+hz), 
Vhat+ VAL Vhraet+ Vhr 2 


eee ap — AE (Wha + Vhz) . 
B=-(Br+B,), d*= at + uz + 6°), 
peers a Oe 
:. AE (Wh vh 
@ _ Ab (Whe + vie) (v2, + urv, +3 +), 


al hip hie alah 
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sF € [2VhrVhi(urve + ure) + (hur + heur)(vr t+ vz)| 
(Vhrhy + Vhthe)(Vhr t+ Viz) 


and 


aC g [2VhrVhL(vRuR 7 ULUL) + (hivR + hrvz)(ur + ur)| 
(Jhrht + Vache) (Wher + Viz) 


5.4 2D Channel Test Problem: Conical Sand Dune 
Test Problem 


To determine which formulation and scheme is the most accurate, we use a conical 
sand dune test problem as discussed by Chesher et al. [2] and De Vriend [6]. The 
test problem consists of a channel of length 1000m x 1000m with dummy initial 


conditions 


Q 


n(e0,0) == 6,90), aa) = a, y,0)’ 


e970) =0 
and the initial bathymetry consists of a conical sand dune, 


- 2 [ x(x—300) - 2 ( m(y—400) j 
B(a,y,0) = ~ (“0”) x (550) if 300 < x < 500, 400 < y < 600, 
0 


otherwise 


and is illustrated in Figure 5.1. The upstream boundary has constant discharge, 


Art -_ hn 


—t,j 0,7? 


Gh =O, (vh)@ yi = (vh)o,; and Be = Boj; 


—1,j 
where Q is a constant and the downstream boundary conditions are 


n+1 n 


n+1 n _ 
Wi and Witji — Wi 


Ww =wrj;, W 


n+l __ 
I+i,j ~~ 


t,-Jj 


where 2,7 = 1 to 5. To obtain the initial conditions, the water flow with the dummy 


initial conditions and boundary conditions are iterated to an equilibrium state whilst 
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Figure 5.1: Initial conditions for 2D Channel Test Problem A (B). 


keeping the bed fixed, where 


lw? — w;.;| < tol 
and tol is the desired tolerance level. The sediment transport fluxes (5.7) are used 
and the porosity is taken as « = 0.4. To ensure the error of the numerical schemes 


do not grow, the variables are non-dimensionalised, 


7 y t h B 

x | i Yy cL is L’ Ee 
~. gl" . vu. , AL . tol . ul 
Pa ge gay ea ee 


where 
L 
L= maxes (larg — tol: lteg = Hol) and T= 7 


denote the non-dimensional coefficients. In this thesis, we take @ = 10 from which 


the initial conditions illustrated in Figure 5.2 may be obtained. Notice that small 
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Figure 5.2: Initial conditions for 2D Channel Test Problem with @ = 10. 


kinks are present in the initial conditions, which do not grow over time. Unless 


stated, the initial conditions illustrated are used for the 2D Channel Test Problem. 


5.4.1 Approximate Solution for Angle of Spread 


De Vriend |6] derived an approximate solution for the angle of spread of the pulse in 
the riverbed. For the conical sand dune test case, the pulse in the riverbed gradually 
changes into a star shaped pattern, which spreads out over time. An approximate 


solution for the angle of spread of the star shaped pattern can be obtained by 
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Figure 5.3: Approximate solution for angle of spread. 


assuming that the riverbed is interacting slowly with the water flow, i.e. A < 0.01, 
and the water flow is in an equilibrium state. This allows the equations of water 
flow, ie. (5.1), (5.2) and (5.3), to be simplified, 


itsizg: OR) SU pw  g(h 4 | eh 
m (huv), + 5 + g(h+ a) =, 


y 
Thus, De Vriend determined an approximate solution for the angle of spread for the 


star shaped pattern, 


tana = OT, — 8T;,’ 
where n 
U | 0g h | 0g 
— —|—1 d T,=-|)—|-1, 
ise) 2 ™ 5 [be 


which is illustrated in Figure 5.3. Here, ¢ denotes the streamwise sediment transport 
flux and U = Vu?+v? the streamwise velocity. For example, if the sediment 


transport flux (3.3) is used in two dimensions, then the streamwise sediment transport 
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flux is 


G= AU" =A (u? +0°)? 
from which we obtain 
i= a | Amo) —-l=m-1 
AU 
and ae 
n [a(AU") 
T, = —m —-l1=-1 
AU Oh 


Hence, the analytical approximation of the angle of spread for the star shaped 


pattern is 


3/3(m — 1) 
tana = ———__,, 
9m — 1 


and if m = 3 then the angle of spread is 


a3 
@= tan! (34) — 21.7867893°. 


5.5 Numerical Results 


In the next few sections, we use the classic 2D Lax-Wendroff scheme, the basic 
scheme (B-FLR) and dimensional splitting scheme (DS2-FLR) with Formulations 
A and C to obtain an approximation for the 2D Channel Test Problem. The 2D 
Channel Test Problem is used with a riverbed that is interacting slowly and quickly 
with the water flow so that Formulations A and C can be compared. Both the basic 
and dimensional splitting schemes use the flux-limited 2D version of Roe’s numerical 
flux (4.14). All numerical schemes in this section are used with Ax = Ay = 20m and 
the minmod flux-limiter. The basic scheme is used with v = 0.4 and the dimensional 
splitting scheme is used with v = 0.8. Formulation A is used with a tolerance level 
of tol = 0.00001 and a cut off time of 15 minutes is used to determine if the water 


flow has reached an equilibrium state. 
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5.5.1 Formulation A: Numerical Results for a Bed which is 


Interacting Slowly with the Water Flow 


We first consider the 2D Channel Test Problem with a riverbed that is interacting 
slowly with the water flow, i.e. A = 0.001. The classic 2D Lax-Wendroff scheme, 
the basic scheme (B-FLR) and the dimensional splitting scheme (DS2-FLR) are 
used with Formulation A to obtain an approximation of the conical sand dune test 
problem at ¢ = 100 hours. The numerical results of the classic 2D Lax-Wendroft 
scheme with Formulation A are not illustrated due to the scheme producing spurious 
oscillations in the results almost immediately, which rendered the scheme unstable. 
From Figure 5.4 and Figure 5.5 we can see that the B-FLR scheme has produced 
smooth numerical results that do not suffer from spurious oscillations. Also, from 
Figure 5.6 and Figure 5.7 we can see that the DS2-FLR scheme has produced smooth 
numerical results that are similar to the B-FLR scheme. However, both schemes 
produced kinks in the numerical results and DS2-FLR scheme has produced a slightly 
different star shaped pattern than the B-FLR scheme. Figure 5.8 and Figure 5.9 
illustrate the approximate angle of spread compared to the numerical results for the 
B-FLR scheme and DS2-FLR scheme respectively. Here we can see that the DS2- 
FLR scheme produced an angle of spread for the star shaped pattern that was closer 
to the approximate angle of spread than the B-FLR scheme. Also, notice that from 
Figure 5.10 we can see that both schemes produced small kinks in the numerical 
results and have increased the total height of the river by a small amount. However, 
the DS2-FLR scheme produced numerical results with fewer kinks present and with 
an angle of spread that was closer to the approximate value. Thus, the DS2-FLR 
scheme is the most accurate scheme. The numerical results of the B-FLR scheme 
and the DS2-FLR scheme are devoid of the spurious oscillations present in the 
numerical results obtained with PISCES, see Chesher et al. [2], in some simulations 


and therefore represent a notable improvement. 
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Figure 5.4: Numerical results for Formulation A using the basic scheme (B-FLR) 
with A = 0.001 and at t = 100h (B). 
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Figure 5.5: Numerical results for Formulation A using the basic scheme (B-FLR) 
with A = 0.001 and at t = 100h (B). 
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Figure 5.6: Numerical results for Formulation A using the dimensional splitting 
scheme (DS2-FLR) with A = 0.001 and at t = 100h (B). 
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Figure 5.7: Numerical results for Formulation A using the dimensional splitting 
scheme (DS2-FLR) with A = 0.001 and at t = 100h (B). 
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Figure 5.8: Illustration of the angle of spread for Formulation A using the basic 
scheme (B-FLR) with A = 0.001 (B). 
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Figure 5.9: Illustration of the angle of spread for Formulation A using the dimen- 
sional splitting scheme (DS2-FLR) with A = 0.001 (B). 
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Figure 5.10: Numerical results for Formulation A using the basic scheme (B-FLR) 
and the dimensional splitting scheme (DS2-FLR) with A = 0.001 and at t = 100h. 


158 


5.5.2 Formulation C: Numerical Results for a Bed which is 


Interacting Slowly with the Water Flow 


Formulation C is now used with the classic 2D Lax-Wendroff scheme, the basic 
scheme (B-FLR) and the dimensional splitting scheme (DS2-FLR) to obtain a 
numerical approximation of the 2D Channel Test Problem for a riverbed that is 
interacting slowly with the water flow, i.e. A = 0.001, at ¢ = 100 hours. From 
Figure 5.11 and Figure 5.12 we can see that the classic 2D Lax-Wendroff scheme has 
produced spurious oscillations in the numerical results. However, from Figure 5.13 
and Figure 5.14 we can see that the B-FLR scheme has produced smooth numerical 
results that do not suffer from spurious oscillations. Also, from Figure 5.15 
and Figure 5.16 we can see that the DS2-FLR scheme has produced smooth 
numerical results that are similar to the B-FLR scheme. Notice that no kinks are 
present in the numerical results obtained with Formulation C, but kinks were present 
in the results obtained with Formulation A. Figure 5.17 and Figure 5.18 illustrate the 
approximate angle of spread compared to the numerical results for the B-FLR 
scheme and DS2-FLR scheme respectively. As with the numerical results obtained 
with Formulation A, here we can see that the DS2-FLR scheme produced an angle 
of spread for the star shaped pattern that was closer to the approximate angle of 
spread than the B-FLR scheme. The angle of spread obtained from using 
Formulation C with the B-FLR scheme and the DS2-FLR scheme are also closer 
to the approximate value than the angle of spread obtained from using Formulation 
A. However, from Figure 5.19 we can see that the numerical results obtained with 
Formulation C have significantly increased the total height of the river and velocity 
from their original values. This is due to all schemes requiring a considerable amount 
of time steps to reach the final computation time of 100 hours and the upstream 
boundary also having constant discharge imposed. Thus, a numerical scheme that 
can be used with large time steps such as an implicit scheme is required to obtain an 
approximation of Formulation C. Also, notice that a few small kinks are present in 
the total height of the river and velocity for the numerical results of Formulation C 
but they do not appear in the riverbed and are not as prominent as with Formulation 


A. As with the numerical results obtained with Formulation A, the numerical results 
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Figure 5.11: Numerical results for Formulation C using the classic 2D Lax-Wendroff 
scheme with A = 0.001 and at t = 100h (B). 
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Figure 5.12: Numerical results for Formulation C using the classic 2D Lax-Wendroff 
scheme with A = 0.001 and at t = 100h (B). 
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Figure 5.13: Numerical results for Formulation C using the basic scheme (B-FLR) 
with A = 0.001 and at t = 100h (B). 
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Figure 5.14: Numerical results for Formulation C using the basic scheme (B-FLR) 
with A = 0.001 and at t = 100h (B). 
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Figure 5.15: Numerical results for Formulation C using the dimensional splitting 
scheme (DS2-FLR) with A = 0.001 and at t = 100h (B). 
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Figure 5.16: Numerical results for Formulation C using the dimensional splitting 
scheme (DS2-FLR) with A = 0.001 and at t = 100h (B). 
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Figure 5.17: Illustration of the angle of spread for Formulation C using the basic 
scheme (B-FLR) with A = 0.001 (B). 
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Figure 5.18: Illustration of the angle of spread for Formulation C using the dimen- 
sional splitting scheme (DS2-FLR) with A = 0.001 (B). 
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Figure 5.19: Numerical results for Formulation C using the basic scheme (B-FLR) 
and the dimensional splitting scheme (DS2-FLR) with A = 0.001 and at t = 100h. 
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of the B-FLR scheme and the DS2-FLR scheme are devoid of the spurious oscillations 
present in the numerical results obtained with PISCES, see Chesher et al. [2], in some 


simulations and therefore represent a notable improvement. 


5.5.3. Formulation A: Numerical Results for a Bed which is 
Interacting Quickly with the Water Flow 


For the second test problem, we use the 2D Channel Test Problem with a riverbed 
that is now interacting quickly with the water flow, i.e. A = 1. The classic 2D Lax- 
Wendroff scheme, the basic scheme (B-FLR) and the dimensional splitting scheme 
(DS2-FLR) are used with Formulation A to obtain an approximation of the conical 
sand dune test problem at t = 500 seconds. The numerical results of the classic 
2D Lax-Wendroff scheme with Formulation A are not illustrated due to the scheme 
producing spurious oscillations in the results almost immediately, which rendered 
the scheme unstable. From Figure 5.20 and Figure 5.21 we can see that the B-FLR 
scheme produced numerical results with no spurious oscillations present. By using 
the DS2-FLR scheme, we obtain the numerical results illustrated in Figure 5.22 and 
Figure 5.23. Here, we can see that the DS2-FLR scheme has produced numerical 
results that differ to the B-FLR scheme. This is due both schemes producing 
completely different star shaped patterns. The DS2-FLR scheme has produced a step 
in the ramp behind the star-shaped pattern whereas the B-FLR scheme produced 
a smooth ramp without a step present. The step present in the numerical results 
of the DS2-FLR scheme was also present when A = 0.001 was used, but the step 
is now larger. The angle of spread of the star shaped patterns also differs. Also, 
notice that both schemes have produced kinks in the numerical results. Figure 5.24 
shows the total height of the river and velocity for the B-FLR scheme and DS2-FLR 
scheme. Notice that both schemes have produced slightly different numerical results 
where the height of the river and velocity has changed from the original values. This 
is due to the formulation imposing a constant discharge on the upstream boundary 
and requiring a considerable amount of time steps to reach the final computation 


time of 500 seconds. 
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Figure 5.20: Numerical results for Formulation A using the basic scheme (B-FLR) 
with A= 1 and at t = 500s (B). 
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Figure 5.21: Numerical results for Formulation A using the basic scheme (B-FLR) 
with A= 1 and at t = 500s (B). 
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Figure 5.22: Numerical results for Formulation A using the dimensional splitting 
scheme (DS2-FLR) with A = 1 and at t = 500s (B). 
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Figure 5.23: Numerical results for Formulation A using the dimensional splitting 
scheme (DS2-FLR) with A = 1 and at t = 500s (B). 
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Figure 5.24: Numerical results for Formulation A using the basic scheme (B-FLR) 
and the dimensional splitting scheme (DS2-FLR) with A = 1 and at t = 500s. 
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5.5.4 Formulation C: Numerical Results for a Bed which is 
Interacting Quickly with the Water Flow 


Formulation A is now used with the classic 2D Lax-Wendroff scheme, the basic 
scheme (B-FLR) and the dimensional splitting scheme (DS2-FLR) to obtain a 
numerical approximation of the 2D Channel Test Problem for a riverbed that is 
interacting quickly with the water flow, i.e. A = 1, at t = 500 seconds. Figure 5.25 
and Figure 5.26 illustrate the numerical results of the classic 2D Lax-Wendroff 
scheme. Here, as with A = 0.001, we can see that the classic 2D Lax-Wendroff 
scheme has produced very poor numerical results due to spurious oscillations 
occurring in the numerical results. From Figure 5.27 and Figure 5.28 we can 
see that the B-FLR scheme produced smooth numerical results with no spurious 
oscillations or kinks present. By using the DS2-FLR scheme, we obtain the numerical 
results illustrated in Figure 5.22 and Figure 5.23. Here, we can see that as with 
Formulation A, the DS2-FLR scheme has produced numerical results that differ 
to the B-FLR scheme due to the star shaped pattern differing. The DS2-FLR 
scheme has also produced a smoother ramp behind the star shaped pattern than 
the B-FLR scheme. The angle of spread of the star shaped pattern also differs for 
the B-FLR scheme and the DS2-FLR scheme. Also, from Figure 5.31 we can see 
that the total height of the river and velocity have remained close to the initial 
values and the results are also smooth with very few kinks present. As with 
A = 0.001, the numerical results obtained with Formulations A and C differed, but 
considerably more with A = 1 due to the sand dune being moved slightly faster with 
Formulation A than Formulation C. This significant difference also occurred in one 
dimension, see Chapter 3, when a large value of A was used and was due to 
Formulation A-CV assuming the water flow is in an equilibrium state, which 


resulted in a constant discharge being imposed. 
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Figure 5.25: Numerical results for Formulation C using the classic 2D Lax-Wendroff 
scheme with A = 1 and at t = 500s (B). 
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Figure 5.26: Numerical results for Formulation C using the classic 2D Lax-Wendroff 
scheme with A = 1 and at t = 500s (B). 
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Figure 5.27: Numerical results for Formulation C using the basic scheme (B-FLR) 
with A= 1 and at t = 500s (B). 
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Figure 5.28: Numerical results for Formulation C using the basic scheme (B-FLR) 
with A= 1 and at t = 500s (B). 
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Figure 5.29: Numerical results for Formulation C using the dimensional splitting 
scheme (DS2-FLR) with A = 1 and at t = 500s (B). 
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Figure 5.30: Numerical results for Formulation C using the dimensional splitting 
scheme (DS2-FLR) with A = 1 and at t = 500s (B). 


172 


Basic scheme Dimensional splitting scheme 


10.05 =" 
Cc 
ne) c 
© g 
® g 
o xo) 
(0) (0) 
8 g 
5 g . 1000 
Paes oe 1000 = 
9.95 : 9.95 + 5 500 
500 as 
500 500 
1000 9 yy jo00 OY 
Xx Xx 
Basic scheme Dimensional splitting scheme 
(04—° 1.04 5 
ry eo 
8 1.02- 8 1.02 - 
ES > 
x x 
3 { Ss 
1000 1 1000 
0 0 : 
500 500 
1000 9 —y 1000 Oy 


Figure 5.31: Numerical Results for Formulation C using the basic scheme (B-FLR) 
and the dimensional splitting scheme (DS2-FLR) with A = 1 and at t = 500s. 
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5.6 Summary 


In this chapter, we extended Formulations A-CV and C to two dimensions and 
adapted the classic 2D Lax-Wendroff scheme, the basic scheme (B-FLR) and the 
dimensional splitting scheme (DS2-FLR) to approximate each formulation. A two 
dimensional conical sand dune test problem was used to determine which scheme 
was the most accurate. The 2D Channel Test Problem was used with a riverbed 
that interacts either quickly, A = 1, or slowly, A = 0.001, with the water flow so that 
the different formulations and numerical schemes could be thoroughly tested. For 
both formulations, the classic 2D Lax-Wendroff scheme produced poor numerical 
results for all values of A and with Formulation A, the scheme became unstable 
almost immediately due to spurious oscillations overpowering the numerical results. 
Formulation A with either the B-FLR scheme or the DS2-FLR scheme produced 
numerical results with no spurious oscillations present, but kinks appeared in the 
numerical results for all values of A. The numerical results obtained with 
Formulation C and either the B-FLR scheme or the DS2-FLR scheme were very 
smooth and did not suffer from spurious oscillations or kinks, which were present 
when Formulation A was used. Notice that the numerical results obtained with 
the B-FLR scheme and the DS2-FLR scheme differ for both formulations. When 
Formulation A is used with the DS2-FLR scheme, a step appears in the ramp 
behind the star shaped pattern, which is more prominent when A = 1, and does not 
appear when Formulation C or the B-FLR scheme is used. Also, the angle of spread 
obtained produced by the DS2-FLR scheme is closer to the approximate value than 
the B-FLR scheme. The numerical results of the riverbed for both formulations were 
similar when A = 0.001 but when A = 1 was used, Formulation A moved the sand 
dune at a faster wave speed than Formulation C. When A = 0.001 was used, the 
schemes with Formulation C dramatically changed the surface elevation and velocity 
whereas Formulation A only changed the values slightly. However, when A = 1 was 
used, the different schemes with Formulation A changed the surface elevation and 
velocity significantly from the initial values whereas Formulation C with the schemes 
produced values very close. This is due to Formulation C with the different schemes 


requiring a considerable amount of time steps to reach the final computation time 
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of t = 100 hours when A = 0.001 but requiring very few time steps when A = 1 
due to the final computation time now being t = 500 seconds. Formulation C with 
the B-FLR scheme uses an average time step of At = 0.72 seconds for all values of 
A < 1 due to the eigenvalues of the water flow being the dominant values. Thus, for 
A = 0.001, the B-FLR scheme requires approximately 500,000 time steps to reach the 
final computation time whereas with A = 1, the scheme only requires approximately 
695 time steps. A constant discharge is imposed at the upstream boundary at each 
time step thus, when a small value of A is used, a large amount of time steps are 
required, which results in the boundary condition changing the surface elevation and 
velocity. Formulation A also produces a different surface elevation and velocity from 
the initial values, but for a large value of A. This is due to the formulation assuming 
the water flow is in an equilibrium state so that the water flow can be approximated 
separately from the water flow and a large morphological time step can then be 
used. For the B-FLR scheme, Formulation A has an overall morphological time step 


that is computed using 
ve vmin (Az, Ay) 


~ 2maxi; (|A" |, |AG])’ 


where 
— 3Agu 


h 


Now, from the initial values we obtain 


a 


(uw? + v?) and \°% = a 


(u? + v’) : 
min;;(\R|) ¥9, max,;(|ul) + 1.08 and max;,;(|v|) + 0.01, 
which implies that 
max; ;(|A"|) © 0.6999A and max;,;(|A°|) = 0.0064805556A. 
Thus, 
yvmin(Az, Ay) 


1.3998A 


and since values of v = 0.8 and Ax = Ay = 20 were also used, we obtain 


At & 


At + 11.4302A7! seconds 
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Hence, when a small value of A is used, the morphological time step is large, i.e. 
A = 0.001 implies that At ~ 11, 430 seconds, but as A > 1, At — 11 seconds. Also, 
Formulation A iterates the water flow to an equilibrium state each time the bed is 
updated, which has a cut off time of t = 15 minutes and results in Formulation 
A requiring a considerable amount of computations when a large value of A is 
used. Thus, a numerical scheme that can be used with large time steps, such as 
implicit schemes, or more accurate boundary conditions are required to obtain a 
better approximation of the different formulations. Hence, Formulation A can only 
be used for small values of A. The DS2-FLR scheme produced considerably more 
accurate numerical results for all formulations than any of the other schemes due to 
the angle of spread of the star shaped pattern and the numerical results obtained 


were very smooth. 
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Chapter 6 


Conclusions and Further Work 


6.1 Conclusion 


Throughout this thesis, we have discussed a variety of numerical techniques for 
approximating the equations governing sediment transport in one and two 
dimensions. In Chapter 1, the equations governing sediment transport were derived 


and two different sediment transport fluxes were discussed. 


Chapter 2 discussed how to obtain an accurate approximation of one dimensional 
systems of conservation laws with source terms. The shallow water equations were 
used to illustrate the different numerical techniques and the accuracy of the schemes. 
The Lax-Friedrichs scheme, classic Lax-Wendroff scheme, MacCormack scheme and 
Roe’s scheme were all adapted to approximate systems of conservation laws with 
source terms. High resolution schemes were also constructed so that the numerical 
scheme satisfies the Total Variational Diminishing property, which ensures no 
spurious oscillations occur in the numerical results. Three test problems were used 
to determine which numerical scheme was the most accurate in one dimension. The 


flux-limited version of Roe’s scheme was the most accurate scheme in one dimension. 


Wars 


In Chapter 3, five different formulations based on either a steady or unsteady 
approach were derived to obtain an approximation of the equations governing 
sediment transport. The flux-limited version of Roe’s scheme and the classic Lax- 
Wendroff scheme were used to obtain a numerical approximation of the different 
formulations. Two test problems were considered and the numerical results compared 
to determine which formulation and scheme was the most accurate. We illustrated 
that the flux-limited version of Roe’s scheme was considerably more accurate than 
the classic Lax-Wendroff scheme, which produced poor numerical results for all 
formulations and became unstable for all variations of Formulation A. Formulation 
A-SF with all schemes produced poor numerical results due to spurious oscillations 
being present. We also showed that Formulation A-NC can only be used for small 
values of A unless the staggered scheme (3.18) was used. We also illustrated 
that Formulation A-CV, which is a steady approach, can only be used for small 
values of A and Q whereas Formulations B and C, which are unsteady approaches, 
can be used for all values of A and Q. Unfortunately, the numerical schemes of 
Formulations A-NC, B and C suffered badly from diffusion due to the schemes 
requiring a considerable amount of time steps to reach the final computation time 


when a small value of A was used. 


Chapter 4 discussed how to obtain an accurate approximation of a system of 
conservation laws with source term in two dimension. The 2D shallow water 
equations were used to illustrate the numerical techniques and the accuracy of the 
different schemes. The classic 2D Lax-Wendroff scheme, the flux-limited 2D Lax- 
Wendroff scheme and the flux-limited version of Roe’s scheme were extended to two 
dimensions. A basic scheme and a dimensional splitting scheme was also discussed. 
A 2D test problem was used to determine which scheme produced the most accurate 
numerical results. The numerical schemes were compared for a 2D test problem and 


the DS2-FLR scheme produced the most accurate numerical results. 


In Chapter 5, the equations governing sediment transport in two dimensions 
were approximated. Formulations A-CV and C were extended from one dimension to 
approximate the equations in two dimensions. The classic 2D Lax-Wendroff scheme, 
the basic scheme (B-FLR) and the dimensional splitting scheme (DS2-FLR) were 
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adapted to approximate Formulations A and C. We illustrated that as with the 
one dimensional case, the classic 2D Lax-Wendroff scheme produced poor numerical 
results for all formulations and became unstable almost immediately with 
Formulation A. Both the B-FLR scheme and the DS2-FLR scheme produced kinks 
in the numerical results with Formulation A but did not with Formulation C, where 
the numerical results were very smooth. The DS2-FLR scheme produced more 
accurate numerical results than the B-FLR scheme due to the angle of spread for 
the star shaped pattern being closer to the approximate value. We also illustrated 
that Formulation A can only be used for small values of A due to the formulation 


being a steady approach whereas Formulation C can be used for all values of A. 


Therefore, in this thesis we have illustrated that formulations based on the steady 
approach can only be used for small values of A and Q whereas formulations based 
on the unsteady approach, with the exception of Formulation A-SF, can be used 
for all values of A and Q. If a large value of A was used, the numerical results of 
the steady approached differed significantly to the unsteady approach. The classic 
Lax-Wendroff scheme produced poor numerical results for all formulations in one 
and two dimensions whereas the adaptations of the flux-limited version of Roe’s 


scheme produced very accurate numerical results. 


6.2 Further Work 


There is a considerable amount of further work that requires investigation. For 
example, Formulations A-CV and A-NC require a more robust and accurate 
approximation of the wave speed of the bed-updating equation. This can be 
considerably difficult to obtain, especially if the sediment transport flux is determined 
by a “black box” approach. In this thesis, all the formulations and numerical schemes 
discussed in one and two dimensions are based on the sediment transport flux 
discussed by Grass [14] and need to be adapted to approximate a general 
sediment transport flux, including the “black box” approach. Numerical schemes 


that can use large time steps, such as implicit schemes need to be investigated with 
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the unsteady approaches so that the numerical results of the unsteady approach 
are less diffusive. Friction, wind resistance and other factors can be included in the 
shallow water equations and the different formulations and numerical schemes 
adapted. For Formulation A-CV, instead of iterating the equations to an 
equilibrium state and then updating the bed, an approximation of the two 
simplified equations (1.3) and (3.2), which were discussed by Cunge et al. [3], 
can be obtained instead. In two dimensions, finite volume schemes and genuinely 
multidimensional upwinding schemes, see Hubbard & Baines [17] and Garcia-Navarro 
et al. [10] for more details, need to be investigated and used to numerically 
approximate the different formulations. With these scheme, instead of the basic 
cartesian mesh being used, different grids can also be used such as triangles and 
even adaptive unstructured grids, where the mesh is determined to conform with 


the geometry. This enables a more accurate numerical approximation to be obtained. 
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